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ABSTRACT 


Two dereverberation techniques are applied to synthetic 
seismic data and their performance in removing water column 
multiples is evaluated and compared. One method is an applica- 
tion of homomorphic deconvolution and the other utilizes 
linear estimation based on a Minimum mean square error criterion. 
The analytical formulations of both methods are discussed. 
Performance is evaluated in terms of three criteria: percent 
of multiple energy removed, percent of signal (reflector) 
distortion, and visual improvement of the data. Results are 
presented which represent the performance of both algorithms 
for a range of environmental and signal processing parameters 
including white noise level, multiple coherence, reflector/ 
multiple overlap, filter parameters and water column travel 
time estimate. The techniques are found to have comparable 
effectiveness on the synthetic data; however, indications are 
that homomorphic dereverberation has greater potential in 
shallow water applications while the linear technique appears 
to be more efficient for deep water data. 
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CRAPIER. I 
INTRODUCTION AND PROBLEM STATEMENT 

The geologic structure of the earth beneath the seafloor 
is most often determined by seismic profiling. The procedure 
generally involves excitation of an impulsive acoustic source 
near the sea surface and recording of the reflected earth 
response with a hydrophone array. Low frequency sound 
penetrates ae bottom and propagates in the substrata with 
reflections occurring at discontinuities in the acoustic 
impedance of the earth. The thickness and density of sub- 
bottom layers may be estimated from the reflected acoustic 
Signal, or seismogram. : 

The earth is modelled as a discrete jerered medium with 
distinct interfaces for most seismic applications. This 
assumed structure, while not strictly accurate, has led to 
good processing results in practical seismic work and has the 
additional advantage of being analytically tractable. We shall 
employ this assumption throughout the present analysis. A 
detailed description of the earth model used iS given in 
Cnavctewme LIT. | 

The amount of energy reflected at a discontinuity 1s 


ideally measured by the reflection coefficient, 


i oe. = al 


ee eC 
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where rr Yo and C11 Cy are the densities and acoustic 


velocities in media 1 and 2. Here the sound propagates from 
medium 1 to medium 2 ana normal incidence has been assumed. 
The assumption of normally incident plane waves is generally 
made in single channel (one receiver) seismology since it 
leads to analytical simplicity and appears to be reasonably 
aeeueate.s on Utilizing this simplified model we. have ignored 
near field effects and spreading losses. These are not of 
major importance in Single channel dereverberation and their 
inclusion would unnecesSarily complicate the earth model. 

It is evident from the reflection coefficient expression 
“Tnees lee reflections occur at points of significant change 
in impedance. The largest impedance discontinuity encountered 
by seismic signals occurs at the water-air interface where R 
is nearly -l. The water-bottom interface is also a strong 
reflector in most cases. Thus, the water column becomes a 
reverberating channel wherein a significant portion of the 
source energy is trapped. Repeated reflections Crom: the bottom, 
or multiples, are received at intervals corresponding to the 
two-way travel time of sound in the water column. Deeper 
reflections from the substrata are masked by multiples when 
their arrivals are nearly coincident in time. Since propaga- 
tion losses are much smaller in water, the energy in the 
multiples is usually large compared to that in the deeper 


reflectors. Thus, water column multiples form an unwanted 





component of the seismogram. 

The reflection coefficient expression indicates that 
all earth layers also introduce reverberation. Except in some 
shallow water situations, internal multiples are generally 
not a serious problem for two major reasons. Most importantly 
the acoustic attenuation in the earth is much greater than 
that in water so that very little internal multiple energy 
is actually returned to the receiver. Secondly, the reflection 
coefficients at earth layer boundaries are usually small 
compared to those at the surface and seafloor so that a 
Bomatively small part sof the aneident energy is actually 
trapped. 

Figure 1 shows a synthetic seismogram with strong 
multiples. Response amplitude is measured on the ordinate 
and travel time in seconds on the abscissa. The first large 
Signal component is the bottom reflection at one second of 
travel time or about 750 meters water depth. The first 
multiple is an attenuated, phase-inverted replica of this 
reflection at two seconds. Note that the return from a 
reflection horizon at two seconds travel time would coincide 
with this multiple and be obscured. The overall periodicity 
of the multiples is apparent in this plot. Actual reflectors 
occur at 1.5, 2.2 and 2.9 seconds. Figure 2 shows the seismic 
environment which would produce such a seismogram. 


The first practical multiple analysis and dereverberation 
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Figure 2 Earth structure which would produce seismogram 
like that shown in Figure l. 
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algorithm was proposed by Backus [1]. His approach was to 
characterize the water peers aS a sharp, ringing filter with 

an impulse response composed of a weighted sum of delayed 
ampulses. The weights and delays are determined by the bottom 
reflection coefficient and the water column travel time 
respectively. This model leads to a three-point operator with 
elements spaced at intervals corresponding to the two-way 

water travel time. Implementation of the Backus filter requires 
estimation of the bottom reflection coefficient and water 

column travel time. Several aspects of the performance of 

this method are discussed in Chapter II. | 

Spatial processing has also been used to reduce 
multiples [2]. Spatial schemes normally require multichannel 
arrays of large physical extent which can be effectively focused 
to discriminate against reverberation. Such systems are widely 
used and quite effective in shallow water but their costs, both 
for hardware and data processing, are very high. Hence, there 
is still a need for time domain multiple removal techniques in 
deep water Situations and for single channel systems. 

Two techniques have recently been applied to seismic 
multiple removal with demonstrated success. The first, an 
inverse filter algorithm based on a tapped delay line model, is 
due to Baggeroer [3], and is referred to hereafter as the TDL 
filter. A tapped delay line is simply a realization of the 


time domain convolution of a signal and a gapped operator [4]. 








Baie 
Results of this method have proven superior to those of three- 
porns Liitering, at least in some deep water situations. 

A nonlinear filtering scheme using homomorphic deconvolu- 
tion has been applied to several aspects of seismic signal 
processing. The use of this method for dereverberation nee 
been demonstrated by Stoffa, Buhl and Bryan[5]. Tne homomorphic 
transformation is essentially a mapping from convolution to 
addition so that, after transforming, deconvolution can be 
accomplished by simple linear filtering. Seismic dereverbera- 
tion appears theoretically to be a very promising application 
of homomorphic deconvolution because of the distinct properties 
of seismic signal components. The method has not, however, 
been fully evaluated or widely used in practice. 

The motivation for this study arises from the disparate 
theoretical mechanisms by which these techniques operate to 
perform the same function. Since analytical comparison is 
not feasible, this functional, comparative approach is 
thought to be the best means of gaining insight into this 
interesting problem. 

The purpose of the analysis is twofold. First, it is 
intended to indicate those factors which have significant 
effects on the performance of each algorithm. The factors 
to be considered are environmental variables and processing 
parameters. These are discussed in Chapter III. Secondly, 


the analysis is intended to point out the relative strengths 





oe 


and weaknesses of the methods by comparison of results on 


Similar data. 


Each algorithm is evaluated for a range of simulated 


processing conditions. Quantitative and qualitative criteria 


are 
the 
fe @ Ye 
ie @ ts 


the 


specified wnich provide a comprehensive description of 
manner in which each signal is affected by processing 
multiple removal. These criteria also ore as a basis 
comparison of results. The scope of the analysis and 


specific performance criteria are discussed thoroughly 


maeemapter IIL. 
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ANALYTICAL FORMULATION AND IMPLEMENTATION 


MoemidlyereaterOormulation of the TDL Filter 

A simple feedback model for the propagation of reverberat- 
ing seismic signals is given by Baggeroer [3]. Multiple 
removal based on this model is then formulated as an inverse 
filtering problem. The dereverberation filter is designed 
using a least squared error criterion and the constraint that 
the filter have a tapped delay line structure. 

The formulation will be developed here from a different 
point of view using Baggeroer's feedback model as a starting 
point. A summary of the feedback model is included for clarity. 

Figure 3 shows the Laplace transform representation of a 
propagating seismic signal. S(s) is the transform of the 
source Signature. The signal first encounters the downward 
travel time delay which corresponds to a phase shift in this 
domain. Hy (s) represents the transfer function of the earth 
beneath the water column including the reflections from layer 
boundaries which comprise the desired information. Internal 
multiples or reverberations between the various earth layers 
are also included in Hy (s)- Another phase shift corresponds 
the return of the reflected signal through the water column. 
P(s) is the feedback gain representing the water column multiple 
mechanism. In most cases this is well approximated by -l, 


which corresponds to the nearly perfect pressure release 
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reflection at the surface. H, (s) includes the observation 
effects such as nonecchone Gandalaen and array tow depth. 
Ambient noise and reciever front end noise are modelled as 
additive white Gaussian noise. 

The overall transfer function is 

—2ST 
R,(s) Hj(s) Hy(s) e ™ 


Be oe 65a (1) 
S(s) 1 - P(s) H, (s) eo 2 Ty 





where at is the one-way water travel time and Ro fs) is the 
received signal without additive noise. 
It 1S apparent that the presence of multiples is due only 
to the denominator of this expression. Thusfar we have 
assumed implicitly that the earth response can be modelled 
accurately as a linear system and that the multiples are 
Saetly Derlodie. The validity of eS ormentone will 
become apparent in the discussion of performance in Chapter IV. 
The obvious task is now to design an inverse filter having 
the form 


-2sT 
F(s) = {1 - P(s) Hy (s) e ] 


Hence, we are required to estimate T and the impulse response, 
hy (t), corresponding to HL (s)- The earth response need not 
be estimated precisely for its entire duration. Estimating 
the dominant energy part of hy (t) is adequate to produce an 


effective dereverberation filter. A typical deep water 
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seismogram has the great majority of its energy concentrated 
Mamie: tirse 200-300 msec of its duration. Effective represent- 
@i10n of this portion of the signal requires about 10-20 
filter coefficients, depending on the bottom and source 
SGharacteristics. 

The transfer function of eauation (1) can be re-written 


im series form as 


Pes) 00 -2nsT 

: =H (s) J (-1)"*" 4, (s) : 
O b 

S(s) n=1 





The received signal then has the time domain representation 


ee mele) ee) Toes emo 


where * represents convolution. This can be rewritten as the 


sum of the primary return and the multiple signal. 


oO 


ry tt) - [s(t x hg (t) | * fh, (t = 27) + Dan, (e-2nn,,) | 
rf) = Dt m(t) 

where 
Dit) <=' s(ojye- hot) = hy (t-2T, ) 


1s the received primary and 
= * =} * ms 
m(t) Se) ha (t) ab, eS) 
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is the received multiple signal. 
The estimation problem, given the feedback model of 
figure 3, is to pe eames DCE tM Chews oresence Of m( =) and 


white noise, w(t). It is convenient to group the unwanted 


elonat Components. 
Mae Ct) te w(t) ; (2) 


Since the unwanted component is an additive one, we 
can consider estimating n(t) and subtracting the result from 
the received signal. We then have the filtering problei 


depicted in figure 4. 


alte) Su ne ears 11 (t) 


me~e £{t) 15 the filter impulse ressense and ce) is the 
minimum mean square error (MMSE) estimate of n(t), given the 
received signal r(t). | 

The number of digital filter coefficients to be estimated 
1s 2TW+l, where T is the effective duration (portion contain- 


ing about 80% or more of the sicnal energv) of hk, (t) and W 


b & 
Peete ee noi andwidthi. Ene -coetfacients will@unen be 


spaced at the Nyquist sampling interval of 1/2W seconds. 
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MicwoOvEImumealoital filter for m(t) wiil have coefficients, 


fue which miniinize 


ip ; 
e=E oe fn (72 = n(i/2u) | (3a) 
i=, | 
wiexre 
+ ¢ | 
Ge 20a) gor [ Gan) = 2, ; (3b) 
kK=1 2W 


The input in (3b) is shifted by the two-way travel time to avoid 

useless filtering of the signal prior to the bottom reflection. 

[i /2w, i ,-/2W] js the time interval over which n(t) is observed. 
SUbSEMEUEING (3%) i2nto (3a) yields 
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fn Rep | Ci-k) 20 - Ro (k/2W+2T_). (4) 


Here we have assumed stationarity over the duration of the 
multiple period. This assumption has led to effective process- 


ing of both real and synthetic data. From (2) 
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Here it is useful to assume (see ref.[3], p.15) that for shifts 
close to 27. Pre Suess Corre lation Of w(t) and r(t) 1s very 
small compared to ew) which will generally have a peak in 
this region. This is equivalent to assuming that the white 
noise has a very short correlation time compared to 2T.- We then 


have 


Ray (KS2Wt2T, ) Rar (K/2Wt2T, ) 


so that (4) becomes 


ats 
) £f, RL((i-k)/2W) = RL (k/2wt+2T. +) (4a) 


Baggeroer has derived equation (4a) by designing the 
Wrenner filter for b(t) with the Peneeeane that the filter have 
a tapped delay line structure, 1.ée. 

2TW 


oc es) fo (tok/2N—20 |) (5) 
k=o . 


When our estimate, nce is subtracted from the unshifted, 
received signal, the resulting overall filter operation has 
exactly eRe Eoume Or (>). This Indirect approach yields ene 
estimator equations of reference [3] without imposing the TDL 
Sernuetunre directly . 

The above derivation also emphasizes the estimator- 
SWotreeeomroGm prea: ction Error structure of this filter.) The 


entire impulse response may be written as follows: - 





= fa 


Diss k 


mars is the prediction error structure for a prediction distance 
of 2 Equation (4a), however, which generates the tf. 
differs from the prediction equations in that the right hand 
Side vector = Rar (tt2T,) rather than Eee Gc a2 oe) As written, 
equation (4a) corresponds to the Wiener filter which produces 
the MMSE estimate of m(t) with r(t-2T, ) as an input. Subtraction 
of this estimate from r(t) whitens only those spectral compon- 
ents which are due to the multiple. 

Peto otineyeeotovedethat the magnitude of the error in b(t) 


meonequal to that in the prediction operator. 


eee ne) 
b(t) = b(t) + (n(t)-n(t)) 
eee) incon (et 


Thusfar the only departures from optimum estimation have 
been the two assumptions of stationarity’ and the relative 


insignificance of R (t+2T ). One further assumption is 


1 
PetCo EOrmactual implementation ef the f£riltern  seiece thar 
Rar (tt2T,,) is a required input which is apparently not measur- 


able from the given data. Baggeroer has observed that for the 


deep water case, which is of primary interest for this method, 
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a (T+2T, )e RLY (t+2T) ) ‘ 


Midte1s, £0r shifts of nearly twice the water travel time, the 
great majority of the crosscorrelation energy is due to m(t). 
Hence, the equations used for data processing are given by 

eg 
, Cela) 7eW) = 2 (i/2ua20 =O; Vee oT 


Having seen the analytical formulation of this inverse 
faltering procedure it is instructive to compare it with the 
Backus three-point method. Processing actual data with both 
filters (ref.[3]) has shown the Backus filter to be significantly 
less effective. Some reasons for this are apparent from the 
foregoing analysis. 

The Backus filter is rigidly dependent on the accuracy of 
two aSsumptions. The first, the assumption of strictly periodic 
multiples, is violated due to the horizontal separation of 
source and receiver. This effect becomes more severe as water 
depth decreases. Since the Backus filter is implemented as 
only three, equidistant operator coefficients it is very sensitive 
tO this lack of periodicity. Even if the statistics of the 
Signal generate very accurate estimates of the bottom reflec- 
tion coefficient the filter structure iS So Simple and rigid 
that proper cancellation will not occur if the multiples are 


Sigmaricantly aperiodic. 
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A second restrictive assumption of the Backus filter is 
that the bottom reflection should be accurately characterized 
Praea single reflection coefficient. All of the statistical 
information available is forced into a single parameter 
estimation scheme. It is apparent that such a filter lacks 
flexibility for dealing with more complicated bottom interaction 
mechanisms. 

The relatively better performance of the TDL filter on 
real data 1S apparently due to its greater inherent flexibility. 
That is, the finite length impulse response, or prediction 
operator, gives the filter a capability for removing reverbera- 
tion effectively in cases where the bottom response is not 
accurately modelled as a weighted impulse. If the bottom has 
a ringing or smearing effect on the incident signal then the 
deconvolution operator must be extended in time. The Backus 
filter, because of its rigid structure, cannot accomodate these 
Situations. The TDL structure provides 2TWtl (usually 10-20) 
parameters which can be varied in the deSign procedure to 
optimize dereverberation of each seismogram. The special case 
of an ideal bottom will generate a filter response which is 
essentially a single spike proportional to the bottom reflection 
coefficient. This result has been confirmed in the analysis of 
synthetic data. In such a case the TDL filter consists 
basically of the first two points of the Backus three-point 


filter. Performance (multiple energy removed) in thesSe cases 
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was found to be essentially independent of operator length. 

The structural flexibility of the TDL filter also gives 
it some potential for dealing with aperiodic multiples. It 
should be noted, however, that the TDL algorithm, like the 
Backus and other classical dereverberation techniques, is 
essentially a correlation-cancellation operation. Consequently, 
increased aperiodicity of multiples can be expected to degrade 


performance in all cases. 


B. Implementation of the TDL Algorithm 

Figure 5 shows a flow diagram of the filter implementation 
used for this analysis. Actual programming was done in 
POmiran lv fon Use On a 32K computer. Re remeing to Ligure ss, 
mmewcOorrelation EUncCtLOnN 1S computed by the standard shift-and- 
add operation with no windowing applied. Results of windowing 
are included in reference [3]. Correlation time is variable 
and may be specified by the operator. The crosscorrelation 
function 1S approximated by the correlation function as discussed 
in the previous section. 

Solution of the filter equations is accomplished by 
conventional matrix inversion. The foeplitz symmetry can be 
exploited for computational savings. Spacing of the operator 
elements ale determined by the estimated signal bandwidth which 
1s specified as an input parameter. 


Actual deconvolution is implemented exactly as shown in 





Figure 5 
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figure 6, 1.e., by means of a tapped delay line or, equivalently, 


convolution with the gapped Goetacon. 


C. Analytical Formulation of Homomorphic Deconvolution 

A homomorphic system is one which obeys a generalized 
principle of superposition and which can be represented as an 
algebraically linear transformation between two vector spaces. 
A detailed description of the theory of homomorphic systems 
is given by Oppenheim and Schafer [6]. This material will not 
be repeated in depth here; rather, we shall discuss the basic 
characteristics of homomorphic systems for convolution with 
emphasis on those properties which facilitate dereverberation. 
Additional discussions of these properties are found in 
references [5], [7] and [8]. 

The usefulness of linear systems for separating additively 
combined signals is due primarily to superposition. Signals 
which are added and happen to be disjoint in the frequency 
domain can be separated by means of an appropriate bandpass 
filter. Homomorphic systems for convolution have a similar 
effect on signals which have been convolved. That is, a 
homomorphic transformation maps the input Signal to a domain 
in which the convolved components may be disjoint. Such a 


transformation is illustrated in figure 7. 
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Figure 7 





== 
The mapping characteristic of this system is intuitively 
apparent. Suppose x(n) is composed of two components which 


have been convolved. 


x(n) = x, (n) * x5 (n) 


The z-transform operation yields a signal with multiplied 


components, X, (2) and X, (2). The logarithm output is 


X(z) = log [X, (2) ] Teh ec [X, (2) ]. 


The inverse z-transform, 


x(n) = x, (n) + x. (n) ; 

preserves the additive combination of the components and yields 
a sequence which is real and stable for a real and stable x(n). 
The sequence x(n) is called the complex cepstrum of x(n). 
Although it is real for real inputs, the "complex" is retained 


to emphasize that it contains both the magnitude and phase 


information from X(z). Hence the complex logarithm is required 
even for real x(n). (We shall omit the modifier here for 
peeVvity . ) 


The cepstrum variable, n, is normally called the quefrency 
(a paraphrase of frequency) or period variable. Filter opera- 
tions in this domain are generally similar to those encountered 
in the frequency domain. Exact subtraction of x, (n) from ay 
followed by computation of the inverse cepstrum (figure 8) 


yields the sequence x(n), exactly, in the time domain. 





x (1) 





It is this property which renders homomorphic processing valu- 
able for deconvolution. 

As an example of complex cepstrum transformation consider 
a Signal composed of a short pulse (2 samples) convolved with 


a decaying, periodic impulse train. 


at° . 


x(n) = s(n) * p(n) 
x(n) = (6(n) + 5 6(nt1))* J (-1)*r*e (n-k 7) 


.K=0 
where 
a aaelerand 1 > 2 


Taking the z-transform, 


CZ) ge (lab) ) (a pe ee 
k=0 


el 


fy 


X(z) = (1 + 2/2)+(1 + Rz 


The logarithm then produces a sum 


nw 


Miz) = toatl. +272) oe tog ties Re 


Ie 





Oa 


Both terms are simply expanded in Laurent series. 
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The z-transforms are easily recognized. 
A n 
s(n) = = ’ Fee ye ae ee 
MZ 
a ee 
p(n) = i ee nok sk = 12. 
k 
x(n) s(n) + p(n) 


an 


We see that s(n) occupies only the negative quefrency 
region and aoa) only the positive quefrency region. Exact 
Geconvolution can be accomplished in this case DY )ZereomncGesrre 
desired half of the cepstrum. 

Mie wecenonre f£Orm Of HNOMOMOrphic SyYSteEmS for COonVvVoLULlenmrs 


eiovwn in figure 9. 


| a = 5 (vin 


Figure 9 
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ieeeechtaracteristic systems D, and Dee areUonOwl eam figures. 7 
and 8 respectively. The system L is a conventional linear 
eyocem. When Lb has a transfer function of 1, then x(n) is 
recovered exactly at the output of pa The choice of Lwaill 
determine the effectiveness of deconvolution for a given invut 
sequence. 

Two further specifications are required to ensure the 
validity and uniqueness of the transformation. The complex 
logarithm is a multivalued function with an infinite number of 


Pbeanches. 


log[X(z)] = oa xe + j(arglx(2) J: ank) for alk ak 


where Arg specifies the principal value of the phase. This 
ambiguity must be resolved while simultaneously satisfying the 
reguirement that ea be a valid z-transform. Note that if 
mn) 1s to be real and stable, X(z) must be conjugate symmetric 
and analytic in an annulus of the z-plane cont@éining the unit 
Gimele. That is, the real and imaginary parts of X(z) must be 
SenmciImuous functions of Zz in the region including the unit 
circle. The imaginary part, arglX(z)]J, can only be made 

Peet nuous bye Unwrabolng. Are lX(2))| dm Suche Way Sebel 
jumps of + 2mk are removed. This unique unwrapping leads toa 
unique, valid X(z) which transforms to a stable, real Cnn 


The requirement of a continuous phase curve poses some computa- 


tional difficulties which will be discussed in the following 





2 
section. 

Regie sd ppatemt from the foregoing description that 
homomorphic systems have the potential for separating convolved 
Signals. One might expect, however, by analogy with linear 
systems -that deconvolution is most effective for signals with 
certain cepstral properties. This is, in Fact, the case and, 
fortunately, seismic signal components are generally amenable 
to deconvolution. Recall that a selsmogram 1s modelled as a 
convolution of a source signature and an impulsive reflector 
series. Reverberation appears aS a minimum phase, periodic 


addition to the reflector series. Thus, we have 


seg ea) eee) ee ma) ( a)! 


where p(n) 1s the source signature, b(n) is the desired signal 
and m(n) iS reverberation. Generalizations can be made concern- 
fiaeene cepstral properties of each component. 

The source signature is, in general, a mixed phase, short 
duration time sequence. It is clear from the definition of the 


Peeraisrorm that any such finite sequence transforms Eo ad 


rational function of z with no poles except at the origin. In 
general 
ee ZF R 
P(z) =Cz H (i-a,z “) WM (1-b.2) ja; |,[b51< i 
i=] j=1 J 


The as and De represent zeros inside and outside the unit cixcile 


respectively. Zz corresponds to a linear phase shift. 





soos 


A m a & 
P (Z) = 16g (C + ) log(l-a,z “by + ) MeN Met 


L=1 jJ=H1 
TOC ee n=0 
A m a.” 
(n) = —— n>0 
; 4 nh : (6) 
g bi.” 
ie n<0 
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Here it has been assumed that the linear phase term is removed 
before computation. ory 1S a two-sided sequence which is 
always of infinite duration but decays faster than l/n. 

Hence, most of the cepstral energy is concentrated near the 
Smerrency Origin. ; 

The reflector sequence is modelled as a train of randomly 
Spaced impulses which may be mixed phase. Stoffa, Buhl and 
Bryan [5] give a general, but very complicated expression 
for the complex cepstrum of such a sequence. Some specific 
examples are given by Schafer [7]. The resulting cepstrum is 
an impulse train with impulses at the time domain impulse 
locations, at all their multiples, and at various other loca- 
tions, both positive and negative on the quefrency axis. Three 
important observations can be made. 

The cepstrum of a minimum phase reflector train contains 


M@meontributlons £Or megative Guefrency. ~Consteéer the soeerat 


Gase Of Equation (6): in which all the 2 are zero. This 
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Borresponds to all zeros being inside the unit circle, i.e., a 
minimum phase sequence. Note that p(n) and s(n) in the example 
are minimum and maximum phase respectively, leading to causal 
and anticausal cepstra. 

Secondly, it happens that no non-zero cepstral contributions 
occur between the origin (first impulse) and the location of 
the second impulse in time, if the time series is minimum phase. 
Therefore, the cepstrum of a minimum phase impulse train, unlike 
that of a general sequence does not have its contributions 
concentrated near the origin. The cepstrum of a minimum phase 
impulse train will always contain a gap equal to that pence 
mer first two time domain contributions. 

Finally, we observe (see Schafer [7]) that a reflector 


series can easily be made minimum phase by exponential weighting. 


r'(n) = w r(n) Iw] < 1 


-k al 
Regist ©. 7 II [1~ (aw) z | f2- (4/7) 2 | 
=] 


1 
The value of w is chosen so that 


ee ok for all Dan 


Weighting of the impulse train is effected by weighting of the 


entire signal since, if 


Sn) = pin) * bin), 
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then 


w s(n) = (w'p(n)}*(w'b(n)). 


Meee that the reflector series may be made minimum phase without 
making the entire signal minimum phase. Very little weighting 
Momnormally required in practice. 

The remaining contribution to the seismic signal is reverbera- 
tion. This component is merely a special case of a minimum 
phase impulse train in which the impulses are periodic. It is 


eCasily verified that if 


ss 
m(n) = } y(n)é6(n-kn_) 
k=] 


then ; 
m(n) = y(n/n.). 


The cepstrum is, therefore, periodic with the same period as 
the reverberation. 

A useful property for deverberation is derived by Stoffa, 
et al.[{5]. The derivation is summarized here because of its 
direct pertinence to seismic processing. 


Consider aq normalized multiple signal, 


m(n) = } (-1)*R*6(n-i27, ) (7) 
si 


where R is the bottom reflection coefficient. Then, as we 
have seen 

(-1) 7 Rt 
wine ag 


RCE 


I 
jm 8 


| 6(n-2iT ). 
oe 
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Subtracting the first} cepstral contributions and transforming, 





, , VG 
m.(n) = min) - Y SER g(n-250 ) 
J . i 
1=1 
M, (oh = y 7 a2 W 
isjtl 7+ 
M.(z) =1 exp [-% 274271] 
J i=jtl i 


Expanding in a power series 
oo co 1 ge ae) Ab 
ne 2 W 


M.(z) = II ) - === 
J isjt+l k=0 kti 


: 23j+k 
= 2 -27 +k © ©O -2T J 
M.(z) = 1+ y (-Rz* tw) oS (-R2 Uw) : 
kel °K eae ee oe 


(8) 

Comparing (7) and (8), the first j time domain multiples have 
been removed completely and the (jt+l)st multiple is reduced by 
i (jt+1). All succeeding multiples are also reduced. Thus, 
removal of only tne first cepstrum multiple would remove the 
first time domain multiple and reduce the second by 1/2, the 
ines Dyes. CEC. . 

Having seen the cepstral properties of each seismic signal 
component the advantages of homomorphic deconvolution are 
| apparent. The source signature and reflector series have their 


cepstral energy concentrated in different regions of the quefrency 
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axis, thus facilitating removal of the source. The cepstral 
Bentri butions of the eee ten occur at predictable locat- 
tions so that multiple removal is possible. 

Due to the nature of the homomorphic transformation, 
jinear filtering of the cepstrum is not the normal convolution 
Operation but a simple zeroing of the unwanted contributions. 
That is, the linear filter is generally frequency invariant 
rather than time or quefrency invariant. The name "quefrency" 
was adopted to reflect this reversal of the customary time and 
frequency filtering roles. 

The foregoing analysis is based on assumptions nie to 
those employed in classical seismic processing. Namely, 
we have assumed that the seismogram consists of a source 
Signature convolved with distinct, impulsive reflectors and 
periodic multiples. It is difficult to predict the sensitivity 
of the overall processing scheme to these assumptions because 
of its complex analytical structure. Hence, various parameters 
have been varied in the performance analysis to obtain an 
empirical measure of this sensitivity. 

Finally, we note that the additive noise was not included 
in the analytical formulation of the homomorphic processing 
scheme. The algorithm is designed to separate convolved 
components and, unfortunately, no effective processing gain is 
achieved over added noise. In practive, additive noise has 


been dealt with through classical bandpass filtering. This 
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performance analysis includes a description of the effects of 
additive noise on homomorphic deconvolution. 
D. Implementation of Homomorphic Deconvolution 

Figures 7, 8 and 9 show the seauence of operations required 
for homomorphic deconvolution. A processing scheme was designed 
EOovimplement this algorithm in Fortran IV en Bie oe Catead. 
computer. A ‘flow diagram of the scheme is shown in figure 10. 
Some aspects of the computation are noteworthy. 

All z-transforms in the algorithm are implemented via FFT, 
Recall from the analytical formulations that z-transforms 
involved a the processing of a real, stable sequence are 
required to have regions of convergence which znclude the unit 
mi~ele. The discrete Fourier transform 1s simply a samoling 
of the z-transform on the unit circle which, for properly band- 
limited signals, is sufficient to specify the Signal completely. 
Data sequences are normally padded with zeros to reduce cepstral 
peiasing, @€.¢d., 2048-point cepstra eemeoneerea fon 1024—peme 
seismograms. Since the cepstrum is always of infinite extent, 

a truncated version always results in some aliasing when 
computation 1S not done recursively. 

The major difficulty in computing an accurate cepstrum is 
fieweompuLatiOn Of a continuous phase Curve. The Kaa  SeQuchices 
ecmeericlii waa led at da race ased On Ene rrequency concenc. 
There is no assurance, however, that this sampling rate is adequate 


tO uniquely specify x(z) = log[X(z)]}. 
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The method used for this analysis is due to Tribolet [9] and 
mo tnought to be quite accurate and efficient. The flow 


G@iagram for this algorithm is shown in figure ll. 
; aX z) oun 

The phase derivative, “ay :«COand principal value, Arg[X(z)], 
are easily computed from the transforms of x(n) and on x(n) 
(See appendix A). The values of the derivative are integrated 
using the trapezoidal rule and the integral output is compared 
with Arg X(z) at each step. If the two values do not agree 
within 


20 Tee at) ode ae 2 eae 


where € 1$ a small positive number, the latest computed value 
© arg X(z), say Aes 1s°discarded.- The io) ite tata matene ia returns 
to the last correct integration value, computes an intermediate 
derivative value, and begins integrating with a step size half 
that of the original grid. The integrate-and-compare process 
is continued at this step size until as has been computed 
correctly or until a comparison fails. Integration is resumed 
at the initial step size in the former case or, in the latter, 
step size is again halved. The number of possible step sizes 
is theoretically unlimited. Tne value of e€ may be adjusted 
by trial-and-error for most efficient integration. The 
adaptive step feature compensates for the undersampling problem 
in a very efficient and accurate manner. 

Linear phase contributions are easily identified and 


removed from the computed continuous phase curve. Having 





wuytzrohte Hutddezmun sseyd AZSTOqQtzAL JO AAeCUD MOTA TT eaznhty 


aSWHd 
GCaddVaMNnn 








(M) d 








oe SNTPA °Ad 
on pue (-m)q 
SReCAHSRAUL aqndwop 












a le 
PQ eTduos f 
N 





-40- 





deqs ptab 
uTeW BUO 
DReCAHSRUL 





>poeyndwoo 
Gm) 





ad 


GOVUOLS 











a 
-*A*qd pue (m)q JO SENTRA aRZeTpeWwrTSQUT A 
Butqyndwoo roxy squtod aatssaoons ore “m--+-Tm 
"S3uroad Plib seule JuUSoeL Pe S2e Fon Fon (*A*d) sentTeA Tedtoutzd 
| N pue (M)qd ‘e@ATRRSATISZC 
a ie Goan a oseug oTdwes-N 3nduf 
eo © |! xo oe @ 
(Taya ‘ (“m)q 





-41- 
computed the log magnitude and phase, the inverse z-transform 
as straightforward. 

Linear filtering is accomplished by zeroing the unwanted 
cepstrum values. One might consider windowing procedures which 
are common in linear filtering, but these were not employed in 
the present analysis. 

The inverse cepstrum computation is completely straight- 
forward since no ambiguities arise in the exponentiation 
process. The final steps are shifting the output sequence 
by the linear phase value and unweighting the shifted sequence, 


1f necessary. 
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DESCRIPTION OF THE PERFORMANCE ANALYSIS 


The purpose of this analysis is to evaluate, both quantita- 
tively and qualitatively, the performance of the TDL and homo- 
morphic dereverberation techniques. Each method is to be 
evaluated for a variety of simulated seismic conditions in order 
to determine those factors which significantly influence perfor- 
mance. The comparative nature of the analysis is intended to 
emphasize the relative strengths and weaknesses of each technique. 

It should be noted that absolute performance figures are 
not Baer ently valuable, especially when obtained from synthetic 
data. The diverse geological and oceanographic conditions 
encountered in marine seismology coupled with the many different 
processing systems currently employed may be expected to yield 
a range of absolute results. The greater value of this analysis 
is to indicate the parameters, environmental and mathematical, 
which can be expected to affect significantly the performance 
of these algorithms. The numbers obtained provide a measure 
of the relative performance of the two methods under similar 
conditions and, in some cases, provide asymptotic performance 
criteria. Synthetic data were chosen so that signal parameters 
eoulad be accurately controlled. 

Three criteria are specified for comprehensive evaluation 


of performance. The first, most direct, measure of effectiveness 





-A3- 

is the percentage of multiple energy removed from the signal. 
This is easily measured by calculating the zero-lag correlation 
of the signal, in time windows spanning only the multiple 
locations, before and after processing. That is, the squared 
amplitude (energy) cf the signal in the multiple region ie 
computed after processing and subtracted from the squared 
emelatude Of the original signal in that region. This 
difference is divided by the original energy of the multiple 
to yield the fraction of energy removed by processing. The 
use of synthetic data facilitates this method of analysis 
since reflectors and multiples can be placed in disjoint regions 
to avoid ambiguity. 

tae Seconda Criterion 1s the amount of signal distortion 
introduced by dereverberation. This is measured by comparing 
reflector energy before and after processing. As before, 
reflectors and multiples must be disjoint for meaningful results. 

Although these two criteria provide an accurate quantitative 
measure of performance they are restricted to situations in 
which reflectors and multiples do not overlap. The overlap 
case is most important in processing real data since the 
Poaitiples then abscure the reflectors most severely. In order 
to judge performance in these situations we must evaluate 
qualitatively the improvement of visual interpretahility. 
This visual enhancement of reflector-to-multiple ratio is our 


faecdsecri terion. Lt iS an imoortant measure in Spite Of 1ts 
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subjective nature since the primary means of seismogram 
analysis 1s visual interpretation. 

The analysis is limited in scope to the single channel 
processing configuration. Although both techniques are 
applicable, in principle, to multichannel processing, the many 
additional variables involved and unavailability of aporopriate 
synthetic data would lead to a complicated extension of this 
analysis. A Single channel treatment is adequate to identify 
the important performance traits of both methods. 

Parameters to be varied fall into two general categories; 
environmental and operational. The environomental parameters 
include noise level, multiple periodicity, multiple-to-signal 
level and multiple/reflector eer eon Or overlap. These are 
wearLed Within ranges which are thought -to be representative of 
ambient conditions normally encountered in marine seismolocy. 
Effects of noise level are considered only for the case of 
white Gaussian noise. The effects of aperiodicity have not 
been evaluated for the TDL algorithm because it is intended 
primarily for deep water use where only the first multiple is 
usually of interest. Operational parameters refer to those 
which can be controlled during processing. These.include filter 
cutoff frequencies, operator lengths, travel time estimates, 
cepstral stopbands, and weichting. Windowing of the correlation 
function is not evaluated although a discussion of this subject 


is contained in reference [3]. 
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The algorithm used in generating simulated earth impulse 
responses is due eOuunerlawtew Ol eA brLert description of 
Theriault's earth model is given here. 

The model is based on the following assumptions: 

(1) An air gun source generates longitudinal pressure 
waves which impinge upon all earth layers as 
normally incident plane waves. 

(eee carth layers are horizontally infinite, 
parallel and homogeneous. 

(3) Abrupt changes in acoustic impedance occur at 
each layer interface and these boundaries are 
characterized by the customary acoustic reflec- 
ion and transmission poet cence 

(4) The earth has a uniform density. 

(5) The water column is a non-attenuating fluid 
with a perfect pressure release interface at 
its surface. 

(6) Each layer is characterized by a transfer 
function, F(w) which represents the attenuation 
ichmetwoveletime delay fOr that tavern. 

These assumptions are incorporated into a lumped parameter 

model. Figure 12 shows a frequency domain model of a two 


layer earth. The FB (w) have tneswtunctional Lorm 


1 5 
exp|-jwx. (9) 


Joa, +1 ce 
} 
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Figure 12 Two layer earth transfer function schematic 
(after Theriault). 
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where Xe Cs and a, are the ith layer thickness, sound speed 


and layer attenuation parameter respectively. These transfer 
Emetions are combined ene a semi-~group property rather than 
the usual frequency domain multiplication. 

ine s-—muttiplier completes a feedback loop around the 
source which generates the water column multiples. Since the 
water column is assumed to have no eee eon the multiples 
appear in the earth response as impulses, and in the resulting 
seismogram as replicas of the source signature reduced by the 
bottom reflection coefficient. 

The above multiple mechanism is inadequate for representing 
effects of incoherent (distorted) multiples and varying bottom 
interaction mechanisms. These effects are introduced by inser- 
tion of water column attenuation which causes the bottom response 
and multiple responses to be of exponential form given by the 
Fourier transform of (9). Thus the bottom response is extended 
in time and each multiple is a distorted version of the previous 
one. Examples of multiple appearance with and without water 
column attenuation are shown in figure 13a and b. 

The topmost looo of figure 12 simulates the effects of 
finite receiver depth t. Any number of layers with the desired 
Parameters can be combined into an overall earth transfer 
function which is easily transformed to yield the earth impulse 
response. 


Synthetic seismograms for this analysis were generated by 
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Figure 13 


(a) Synthetic seismogram with coherent multiples. 
(Sy Synthetic Setsmogram with dicverteq miitipice 
due to water column attenuation. 
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convolution of various Theriault impulse responses with air 
gun signatures obtained from actual at-sea MeecOmaiigccs) 
typical signature is shown in figure 14. 

It is important to note that these seismograms contain 
all water column multiples and internal multiples of all 
orders. For realistic parameter values earth attenuation 
characteristics usually render internal multiples negligible. 

Finally, we note one drawback of uSing this model for 
multiple removal analysis. The algorithm produces multiples 
which are exactly periodic. This periodicity gives these 
seismograms strong correlation characteristics which are not 
mMewally encountered in practice. The lack of periodicity in 
actual seismograms is due to the horizontal separation of the 
source and receiving array. The difference in travel paths 
arising from this separation is illustratéd below for a primary 


meturn and first multiple. 
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Figure 14 Air gun signature (obtained by at-sea recordings) 
convolved with synthetic earth response functions 
to form seismograms. 





Se = 
For a known separation and water depth the travel time 
eitference can be easily eanlenratean Tie Sswetieet becomes 
minimal in deep water where cos a and cos 8 are approximately 


eeuat tO One. 
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CHAPTER IV 


RESULTS 


peeeeentroduction 

The results of applying the multiple removal methods 
described in Chapter II to synthetic data are described here in 
terms of the criteria of Chapter III. Performance based on the 
First two criteria is emphasized because it can be quantified 
Bren mOre accurately. Specifically, the reflectors and 
multiples have been positioned at distinct locations in most 
cases so that the amount of enercsy removed from reflectors and 
multiples can be measured without ambiguity. This emphasis 
leads, however, to a certain lack of realism in several of the 
Semeete Gata plots. Separation of this kind in an actual 
seismogram would, of course, eliminate the necessity for 
multiple removal. Hence, several cases of interfering reflectors 
and multiples are also shown. Althouch these are not amenable 
to quantitative analysis they can be judged on the basis of 
ehewtnird Criterion, viz., improvement of visual record quality. 

A second deviation from normal processing conditions has 
been ee i ied £O compare effectively the performance of the 
two methods. Homomorphic dereverberation is accomplished in 
meactice (see [5])} concurrently with source deconvolution, 
where practical, since both operations are simply performed 
after the cepstrum has been computed. This leads to a 


considerable amount of energy removal at each multiple and 
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reflector location which is due to source deconvolution alone. 
Although this may fede seo IMprOvenenttmonmtme Hecord, the crrtenien 
of energy removed does not accurately measure dereverberation 
performance in the same way it does for the TDL results. For 
example, source deconvolution might typically lead to 903% 
reduction of multiples and 90% removal of reflector energy. 
This gross change in configuration cannot be effectively 
compared with TDL processing on the basis of energy removed. 
enice,, Source removal has not been accomplished in most cases 
tested. The results of this "partial" processing can be 
quantitatively evaluated and easily compared with TDL results. 
Cepstral filtering which includes source removal, thus retain- 
ing only high quefrency energy, Ae referred to as "longpass 
filtering". Several examoles of longpass filtering are included 
and interpreted in terms of the third criterion. 

The performance of each method is discussed individually 
for various processing situations. A direct comparison of the 
performance of both methods for the same data is also included. 


The comparison is extended in Chapter V. 


Bs Results of TDL Dereverberation Performance 


i. Opveretor Length, Muitiple Distortion and Multiple-to- 
Siena Ratio 


Bicentnboer Of tao gains required for optimum multipole 
removal was found to be highly signal-dependent. Recali from 


Chapter II that the tapped delay line is essentially an 
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estimate of the high energy portion of h(t), the earth impulse 
response. in most applications the reflected bottom response 
1S perenant. The Moor oncwieen DrOdmeceE OL Es) Fesponse 
would then be expected to govern the filter length requirements. 
Measured results confirm this. 

Seismograms containing exactly impulsive (one sample) 
bottom responses exhibited little or no variation of performance 
with filter encrene Typical filter impulse responses for a 
seismogram of this type are shown in figure 15. The reflection 
coefficient in this case is 0.3. The first tap gain is close to 
-0.3 in each response, as would be expected from the Backus 
formulation in which the second Operator point is an estimate 
Sieihe bottom reflection coefficient. Increasing filter length 
can be seen to cause variation in the "estimation" of the reflec- 
tion coefficient. Figure 16 shows energy removed vs. filter 
length for this seismogram. Multiple energy removed decreases 
slightly with increasing filter length. Reflector distortion 
is nearly constant at low values (6% for one and -2% for the 
Benet). The signal used in figure 16 is shown in figure l/a. 
The processed result shown was obtained using only one tap. 

Introduction of a non-impulsive multiple mechanism was 
found to produce a marked dependence of performance on operator 
length. Synthetic seismograms were generated with finite 
length bottom responses, resulting in distorted multiples 


which are not simply weighted replicas of the source signature. 





| -55- 


| 


(a) 


(b) 


Ce) 


Figure 15 TDL impulse responses for three operator lengths. 
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Figure 16 Operator length vs. performance for a signal with 
with impulsive multiples. 
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Figure 17 (a) Seismogram with coherent multiples. 
(b) Result of TDL processing with a one point 
Gecrarsou. 
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Figure 18 shows filter performance vs. operator length 
for three seismograms having different degrees of distortion 
in their bottom reflection mechanisms. This distortion is 
equivalent to extension of the bottom impulse response in time. 
mie Signals associated with curves (1), (2) and (3) are shown 
mmeergures 19a, 195 and 19c respectively. The increasing 
distortion of the multiple at 2.0 seconds is evident. Curve (1) 
corresponds to a signal with very slight multiple distortion. 
Operator length is seen to have no effect on performance. 

The seismogram corresponding to curve (2) has a bottom impulse 
mesponse which 1S Significant for T = .03 seconds. The tap 
spacing in this case is 1/2W = 4.88 msec, so that 2TW = 6.1, 
and six or seven tap gains should be adequate if the signal 
has been properly sampled. Reference to curve (2) confirms 
that increasing the filter length beyond seven does not improve 
performance. Filters of fewer than seven elements yield 
monotonically decreasing performance. The bottom response 
mereeurve (3) iS Significant for .045 seconds SO Ehat, £Omsene 
same tap spacing, 2TW = 9.2 and we anticipate that nine or ten 
taps will he adequate. This is, in fact, the case. 

Peete 2Oushows the 5, 9 and 15 point filter ampulse 
responses associated with curve (3). There is an observable 
convergence to a fe = shape which 15 tne@ecrualmiuncetonal 
form of the synthetic bottom response. Figures 20a and 20b 


exhibit a "diverging tail" effect which was found to be common 
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Figure 18 Performance vs. operator length for three signals 
VWEeMeDOmcOom Interaction Elmes varying from.) 
impulsive to (3) .045 seconds. 
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Seismograms associated with the curves of figure 
18. Caeeeumave C1) s (pyrreeurve. (2)2 (Cc) Veurve= tae 
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Figure 20 Filter immulse responses associated with Curve (3) 
Orerrcunre 18. (a) 5 pointsmt)) 9 pont; 
(ej > DOInks. 
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men the specified filter length is too short. This type of 
phenomenon occurs in some numerical approximation methods when 
an inadequate number of terms is specified. 

fimo apparent from the above results that the Gffect of 
Operator length is closely related to the bottom interaction 
mechanism. As discussed in Chapter II, the filter should be 
an estimated replica of the bottom impulse response when the 
interaction process can be accurately modelled as a convolution 
of the source signature with the bottom response. Figures 
15 and 20 are good examples of this behavior. 

In the cases summarized in figure 18 performance increases 
as multiple distortion ee eeeses This need not be true in 
ceneral since bottom interactions may become very complex. 


mne Slowly varying Poe: 


responses lead to operators which 

Meweneal greater Cancellation effect as they are extended. A 

higher bandwidth bottom impulse response might not exhibit 

this behavior. As it was not possible to include more complicated 

bottom responses in the earth model used, these effects were 

not investigated further. | 
Reflector distortion was found to be nearly constant for 

pede ritdter lengths tested. The reflector at 2.7 seconds was 

Sscentially undistorted in all cases. The 3.5 second reflecror 


Mac an averege distortion of 7%. Figures 21-24 show some 


Processed results for the signals in figure 19. Hach of the 
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Figure 21 (a) Seismogram with very coherent multiples. 
(b) Result of TDL processing with three taps. 
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Figure 22 
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(a) Seismogram with moderately distorted multiples; 
bottom response is Significant wtore05> second. 


(b) Result of processing with 7 taps. 
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Pigure 23 (a) Seismogram with considerable multiple distortion; 
bottom response is Significant for .045 seconds. 
(a Result Of TDL Drocessing with 11 veaps-. 
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first three figures (21, 22 and 23) show signals before and 
after processing with the fewest number of tans required to 
achieve optimum performance. Visually, multiple removal is 
mmccst Complete in each case. The signal of figure 23a is 
Meee in figure 24 after processing with 3, 5 and 7 taps. 
The improvement from figure 24a to 24c is obvious. 

Multiple-to-signal ratio (defined here as the ratio of 
energy in the first multiple to that in the largest sub-bottom 
reflector, abbreviated MSR) was found to be of little importance 
in most cases of interest. For Signals in which multiples are 
large enough to be a problem (comparable to, or larger than 
smaller reflectors), the multiple dominate the crosscorrelation 
function so that an effective filter is generated. For these 
cases the performance was found to be insensitive to the width 
of the time window used for correlation. In the relatively 
less interesting case of signals with small multiples, the 
performance is greatly Gependent on the choice of correlation 
window. If large reflectors are included in the window, perfor- 
fiamee 1S adversely affected because of the large reflector 
contribution to the statistics. For some cases of interest 
this effect may be a consideration in chooSing an appropriate 
correlation window. Inclusion of ieee reflectors should be 


avoided. 
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Epicure 24 


Results of processing seismogram of figure 23a with 
inadequate TDL lengths (a) 3 taps (b) 5 taps (c) 7 taps. 
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2. Water Travel Time Estimate 

It was noted in Chapter II that an estimate of the water 
column eee time is required for implementation of the TDL 
algorithm. This is accomplished in practice by various methods 
including visual estimation, energy detection and correlation 
techniques. The estimate appears in the filter design equations 
as the minimum shift of the Evoseeo Econ tion. cunetons Rar’ 
This estimate may also be identified as the prediction distance 
when the operator is interpreted as a prediction error filter. 

The actual estimation of water column travel time was not 
investigated in this analysis. The effects of travel time 
Beeination on filter performance were, however, considered. 

Figure 25 shows the results of water travel time estimation 
errors for the ideal case of a signal with impulsive multiples. 
The unprocessed seismogram, shown in figure 26a, contains 
Perilectors at 2.7 and 3.5 seconds and a strong multiple at 2.0 
seconds. The actual two-way travel time in this case is 1.0 
second. The strong Sa Madari ty petHeCTmene bottom reflection 
and multiple is apparent. 

First multiple energy removed is very sensitive to travel 
time estimation, with an error of + 10 msec resulting in a 
performance degradation of about 20%. This effect is analogous 
to that observed in matched filter receivers in that the 


coherent signal components exhibit very high correlation for a 


small range of lags. In this case the strong coherence and 
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Figure 25 Effects of water column travel time estimate on 
multiple removed for a signal with coherent multiples. 
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Figure 26 


(a) 
(b) 


(a) 


—— Yn oe 


(b) 





Unvrocessed seismogram with impulsive multinles. 
TDI, result based on Correct Water “Erave. time 
estimation. 





-~7l- 


Rall R2 


With $n nr ann 


(c) 


| 


Wier 


(d) 








Figure 26 (cont'd) -(c) Result of TDL processing based on an 
early water tieaveol Emmewese abo 
(c) Same processing based on a late estimate. 
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multiple periodicity yield the sharp correlation peak for lags 
near the water column ee time, which is also the multiple 
period, 

The second and third multiple performance peaks occur 
faet.08 and 1.04 seconds respectively, and have values of 243% 
and 60%. These shifted and reduced peaks are caused by 
megttization effects. For such an extremely coherent signal, 

a multiple position deviation of one sample (due to sampling 
interval round-off) can cause the cancellation operation 

mowbe severely affected. In this case the second and third 
multiple performance maxima are due to secondary aneescoteeiaten 
peaks introduced by the periodic oscillations in the source 
Signature. These peaks occur at intervals approximately 

eaval to the period of the hasic source (alr gun) frequency. 

The three performance peak values correspond closely to the first 
three air gun pulse amplitudes. 

Figures 26 b, c and d show the results of on-time, early 
and late travel time estimates respectively. In the first case 
the multiple has been effectively removed while the early and 
mane estimates lead to multiples which are still significant 
after processing. 

The effects of travel time estimation error on the 
operator impulse response are seen in figure 27. The optimum 
estimate yields an operator with a large peak at the origin, 


nearly equal to the bottom reflection coefficient (.2), anda 
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Figure 27 TDL filter imculse responses based on travel time 
estimates which are (a) 40 msec early, (b) correct, 
and (c) 20 msec late. 
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smaller peak at the source pulse oscillation interval. The 
initial spike is due to the largest peak in the crosscorrelation 
function (approximated Here by the correlation function) which 
occurs at a shift equal to the two-way travel time. The smaller 
Beak in the filter is introduced by an additional shift of one 
source period. The early (by 40 msec) estimate still allows 
observation of the crosscorrelation maximum and yields the 
mareted Spike of figure 27a. It is evident in figure 25 that 
performance is reasonably good for estimation errors less than 
about .07 seconds, which is the length of the tapped delay 
line. For greater travel time errors the range of correlation 
lags’ comptted does not include the peak; hence, the filter is 
Beerreettve. Late travel time always results in poor perfor- 
mance since the peak is not observed. Figure 27c shows a 
typical impulse response due to late travel time estimation. 

The seismogram of figure 28a contains a small reflector 
at 2.2 seconds and a larger one at 3.0 seconds. The first 
Mmoltipole partially overlaps the smaller reflector and the second 
multiple coincides with the larger. Figures 28 b, ¢ ane 
demonstrate how travel time estimation errors can lead to 
ambiguous results. Figure 28b was processed using the correct 
mravel time, resulting in good resolution of both reflectors. 
The early and late estimates lead to the signals of figures 
28c and d in which the smaller reflector cannot be clearly 


resolved. 





Figure 28 
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Seismogram with overlapping reflectors and 


multiples. 
Result of TDL processing based on a correct 


travel time estimate. 
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Result of TDL processing based on a 
40 msec early travel time estimate. 
(d) Same processing based ona20 msec late 
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Less coherent seismograms exhibit less sensitivity to 
travel time estimation as shown in figure 29. Curves (1), 
(2) and (3) were generated using the seismograms of figure 19 a, 
b and c respectively, which contain increasing distortion in 
their multiples. The MSR is considerably lower than in 
figure 25 so that the later multiples are not clearly visible. 
The crosscorrelation peaks are considerably broader than in 
the impulsive multiple signals. Curve (3), the least coherent 
Signal, exhibits the lowest sensitivity to travel time estimates. 
Increasingly sharper peaks are evident for curves (2) and (l). 
The performance degradation due to over-estimation is precipitous 
in all cases although the rate of fall-off is related to 
Signal coherence. 

Reflector distortion (not shown) was negligible in all 
cases for the first reflector and nearly constant at 7% for 
the second. 

Examples of processed signals for curves (2) and (3) 


are shown in figures 30 and 3l. 


Beerialtrtple Periodicity 

The effects of aperiodic multiples on TDL dereverberation 
performance were not investigated in this analysis since the 
algorithm is designed primarily for deep water use, where only 
one Mbt ole 1s significant in most applications | Successful 


employment of this technique in shallow water is limited since 
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Figure 29 Effects of water column travel time on three signals 
with varying degrees of multiple distortion. 
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Figure 30 Processed scismoqrams associated with ficure 1.9b 
and curve (2) of figure 29. (a) Travel time 
estimate 40 msec early (b) 40 msec late. 
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Figure 31 Processed seismograms associated with figure 19c and 
curve (3) of ficure 29. (a) Travel time estimate 40 


msec early (b) 40 msec late. 
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Ene approximation of the crosscorrelation, Rar! by the 
merrelation, ee near the multiple onset is not generally 
valid. Reflector energy in the multiple regions is usually 
Significant in shallow water seismograms. 

It is apparent from the TDL formulation and the discussions 
of the preceding two sections that the TDL algorithm has 
some capability to compensate for aperiodicity in which the 
deviation is of the order of the tapped delay line length. 
For greater aperiodicities the TDL filter will not be effective 
on later multiples. Figure 25 indicates that even slight 
deviations can be very detrimental in later pe teate removal. 
Dynamic corrections can be applied as suggested by Backus [1], 


however, these are beyond the scope of this study. 


4. Peflector/Multiple Overlap 

Figures 32-36 show seismograms before and after processing 
for several cases in which multiples and reflectors overlap. 
A brief.description of each situation is given here. 

The seismogram of figure 32a has reflectors at 1.5, 2.1 
and 3.0 seconds, and distorted multiples. Significant energy 
near 2.0 seconds is removed by processing but a clearly visible 
response is still present. Other regions of the signal are 
not visibly affected. The visual resolution is not significantly 
improyed in this case, however, stacking of successive shots 


after multiple removal could be employed to enhance the reflector 
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Seismooram with reflectors at 1.5, 2.1 and 3.0 
seconds (a) before and (b) after TDL processing. 
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in this situation. The success of stacking would require 
that the reflector be more coherent than the residual multiple 
after dereverberation. 

tr figure 33 reflectors occur at 1.5, 2.02 and 3.0 seconds. 
Reflectors may be expected to have a strong influence on the 
crosscorrelation function in this situation of extreme overlap. 
The result 1s nearly complete removal of the second reflector 
with the multiple. The last reflector is clearly visible at 
.3.0 seconds. In this example the tapped delay line length 
(.054 second) is greater than the reflector/multiple separation 
(.02 second). The bottom impulse response is significant for 
-035 second so that it cannot be estimated accurately with an 
Operator which is shorter than the reflector/multiple separation. 
Even if the bottom resvonse were much shorter, the effect of 
the large reflector in the crosscorrelation window would lead 
to degraded performance. This example shows worse degradation 
than would be expected in practice because the reflectors are 
very ccherent in this synthetic data. The less coherent 
reflections in deep water signals would generally be less 
susceptable to removal. When reflector energy is Sree 
in the crosscorrelation, however, and TDL length is greater than 
reflector/multiple separation, the filter design algorithm 
essentially interprets the reflector as part of the multiple 
and tries to remove it. 


Figure 34 contains reflectors at 1.6, 1.95 and 3.0 seconds. 
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(a) 


migure 33 Seismogram with reflectors at 1.5, 2.02 and 3.0 
seconds (a) before and (b) after TDL processing. 





Figure 34 
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(a) 


Seismogram with reflectors at 1.6, 1.95 and 3.0 
seconds (a) before and (b) after TDL processing. 
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These locations @re slightly different from those of figure 33, 
Significant energy 1s present at the second reflector location 
after processing.j Since the reflector precedeS the first 
multiple in this case it has a less significant effect on the 
crosscorrelation function than the 2.02 second reflector in the 
previous example. The variability of performance on individual 
seismograms (of similar structure) with overlap again suggests 
the use of stacking to enhance reflectors. 

The reflectors of figure 35 are situated as in figure 34; 
however, the first reflector and multiple are significantly larger 
in this case. This seismogram is more typical of a shallow 
Water return. The processing is quite effective on the first 
and second multiplies. The aperiodicity of reverberation in 
most shallow water seismograms would lead to much poorer 
performance on later multiples, in practice. 

Figure 36 illustrates a situation where the overlap is 
not severe but visual resolution is impaired by the first and 
Pecema multiples. Reflectors occur at 1.55, 2.4 and 3.25 
seconds. Processing results in effective reduction of both 
multiples and sigmificantly improved resolution in the second 
and third reflectors. ; 

We conclude from these examples of reflector/multiple 
overlap that visual improvement due to dereverberation varies 
widely, depending upon several aspects of the individual signal 


structure. The bottom impulse response, the extent of reflector/ 





Figure 35 
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(a) 


Seismogram with reflectors at 1.6, 1.95 and 3.0 
seconds (a) before and (b) after TDL processing. 
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Figure 36 Seismogram with reflectors at 1.55, 2.4 and 3.25 


seconds (a) before and (b) after TDL processing. 
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multiple separation, the relative coherence of reflector and 
multiples, andm eres relative sizes of signal components all 
appear to affect performance. Each of these factors influences 
Ene erfectiveness of the crosscorrelation operation in identify- 
ing energy which is specifically due to the multiples. 

The visual improvement due to dereverberation can best be 
determined by analysis of continuous seismic profiles (presenta- 
tions showing mamy seismograms side-by-side) since coherent 
‘residual energy, if present, becomes apparent as visual trends 
in the data. The effects of stacking can also be observed in 
this format. It was not practical to generate such a profile 
with synthetic data but the single-shot results indicate that 
TDL dereverberation is Be ean useful for enhancing the 
visibility of reflectors which are partially masked by 


multiples. 


5. Additive White Noise and Filtering 

White Gaussian noise was wenetnee in the following manner. 
First, a set of umiformly distributed random numbers was 
obtained using a standard digital routine. A set of approximately 
Gaussian numbers was then formed by summing separate groups of 
twelve of the uniformly distributed Peron The resulting set 
was weighted to obtain the desired standard deviation. Figure 
37 shows a seismogram with two different levels of added noise, 


each lowpass filtered at 50 Hz. 
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Figure 37 Seismograms with added white noise, lowpass at 50 Hz. 
(a) g=50 Ge} g=100. 
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Figure 38 illustrates the performance of the TDL filter 
for various noise levels on a seismogram having disjoint 
reflectors and Penta lest The maximum noise standard deviation 
shown corresponds to a SNR of 0.7 dB. No bandpass filtering 
was done prior to dereverberation. The results presented here 
are typical of the cases tested. Multiple energy removed 
decreases monotonically with increasing noise level; almost 
linearly in the mace of the first multiple. The later multiples 
are more severely affected. As the noise level becomes comparable 
to their amplitudes, the filter becomes ineffective in removing 
them. This behavior is due to the effect of increasing incoherent 
energy in the Signal. The correlation-cancellation operation 
is designed to detect coherent, periodic components. These 
become increasingly masked as more noise is added. For this 
reason seismograms with larger multiples exhibit less sensitivity 
to noise. Seismograms of similar structure (1.e., exactly the 
same reflector and multiple locations) whose multiple energies 
differed by a factor of twenty were Rend to exhibit consider- 
able dereverberation performance differences in high noise. 
Processing of the signal with the larger multiples resulted in 
25% greater removal of first multiple energy. 

Lowpass filtering before dereverberation can lead to a 
Significant improvement in TDL filter performance in some 


Signals. The plotted results of lowpass filtering two noisy 


Signals at various frequencies are shown in figure 39. 
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Figure 38 Effects of noise level on TDL performance for a 
signal with coherent periodic multiples. 
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Figure 39 TDL performance vs. lowpass filter cutoff frequency 
for signals with coherent (lower curve) and distorted 
(upper curve) multiples. 
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(A third order Butterworth filter was used throughout. The 
visual difference in two filtered signals is shown in figure 40. 
The Peet is lowpass eitered at 90 Hz and the second at 30 Hz.) 
The seismogram used to generate the upper curve contains 
Multiples which are considerably distorted while the lower 
curve is due to a more coherent signal. Decreasing the pass- 
band from 90 to 30 Hz produces a 123% eee ence Der rormlance 
in the second case but the less coherent Signal is relatively 
insensitive to filtering. Several other signals evaluated 
were found to exhibit the same behavior. The phenomenon is 
apparently due to the averaging, or smoothing effect of the 
operator in the distortion case. Recall that signals of this 
kind tend to have more extended waveforms in their TDL operators 
foeertigure 20). Convolution of such functions with a noisy 
Signal "smoothes out" the noise because of the significant 
Operator length and the random fluctuation of the noise. This 
has the twofold result of better overall performance in noise 
and less sensitivity to removal of higher frequency noise enercsy. 
Performance degenerates for filter cutoff frequencies below 
30 Hz because the spectral energy of the.signal itself is 
concentrated in this range. Figure 41 shows the distribution 
of energy in the frequency domain for the seismogram associated 
with the upper curve of figure 39. Most of the energy is 
concentrated between 5 and 30 Hz. 


Figures 42-45 illustrate the visual improvement of noisy, 
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Figure 40 Seismograms with equal white noise levels, lowpass 
frlbteredrae (aj 00) He vaiic(o) mo Oma 
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Figure 41 Squared magnitude of the Fourier transform for the 
Signal associated with the upper curve in figure 39. 
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Figure 42 (a) Unprocessed seismogram with added noise, F=50 ize 
(b) Result of TDL processing. 
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Figure 43 (a) Seismogram of figure 42a with high noise, F=50 Hz 
(b) Result of TDL processing. 
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(a) Seismogram with low noise, PF =50 HZ. 
(b) Result of TDL processing. 
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Figure 45 (a) Seismogram of figure 44 with high noise, F=50 Hee 
(b) Result of TDL processing. 
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Figure 46 Operator impulse responses for three different 
noise levels on Same Signal with coherent multiples. 
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lowpass filtered (50Hz) signals after processing. Two different 
seismograms are shown, each at two different noise levels, 
before and after processing. It can be seen from these figures 
that no serious reflector distortion occurs due to noise or 
filtering. 

The presence of white noise has afpreatetenie effect on 
the appearance of the TDL operator. The impulse response itself 
becomes noisy due to the added incoherent energy and a bias is 
introduced into the large cancellation peak. An example of this 
is shown in figure 46 for a signal with coherent multiples. 
The variation in the first tap gain is significant and the change 


in appearance of the overall impulse response is appearent. 


C. Results of Homomorphic Dereverberation Performance 

iemilntroguction 

The theory of homomorphic dereverberation is based on the 
properties of periodic minimum phase impulse train cepstra. The 
incorporation of more realistic effects, such as aperiodicity, 
distorted multiples, and non-impulsive reflectors, leads to 
analytical intractability in most cases. Hence, in presenting 
ehe Beperimen tal results, we view the more complex signal 
structures in terms of their relation to the ideal signal models 
of Chapter II. In so doing we hope to proyide a basic reference 
for interpreting observed phenomena, and to provide a stronger 
intuitive picture of the mechanisms at work in homomorphic 


dereverberation. 
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All cepstra used in this analysis are 2048 points in 
length. Seismograms contain 1024 points except in cases where 
resampling was required. Zeros have been appended as necessary. 
The 2048 point length has been used throughout the analysis for 
Uniformity and reduction of aliasing. In most practical seismic 


processing, shorter cepstrum lengths are adequate. 


2. Multiple Periodicity 

It was shown in Chapter II that later components of a 
periodic, minimum phase impulse train can be significantly 
reduced by removal of the first one or two cepstral contributions. 
This property is potentially valuable for dereverberation, 
especially in shallow water cases where several strong multiples 
may appear in the data. In order for this property to be useful 
it must be relatively insensitive to (at least) slight 
aperiodicity since actual reverberation is not exactly periodic. 
The result of Stoffa, et al. summarized in Chapter II is 
extended here to rapidly decaying, aperiodic, minimum phase 
impulse trains. 

Consider removing the first j cepstral Foner butions of 
ara, the cepstrum of a minimum phase impulse train with 
contributions (-R)2 at n., i 1.045 eee and (Ri < 1s 


We then have 
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Transforming, 


Exponentiating, 


co Q 
M. (z) = exp } ee i 


| - , 
M. (z) II exp NGUE o TR 
J R=jtl g 


This may be expressed as a power series. 
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J R=j+1 i=O ile 
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(10) 
The rapidly decaying coefficients of gen combine multiplica- 
tively to yield the time domain impulse coefficients of oe 
iPmeumy the first cepstral impulse 1s removed (j=1), the 


largest value attained by the third coefficient in brackets is 


(-R) 22 , 


22 ace 
R=2 8 
which is insignificant for reflection coefficients of interest. 


All succeeding terms in (10) will be vanishingly small. We 
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then have the approximate ae 


M.(z) =I f i ca 2 2s! 
7 Q=jt] |. 


Expanding this product we obtain 








RS Gy Mee 00 (2j+k) <n... 
Mi(z)=1+j) [® 2 9K + Fc (eR) Zin 
J k=] J 3 
which transforms to 
| Se ee Jt2 
m.(n) = é&(n) + (~R) « 6é6(n-n +1) + (~R) « §6(n=n ggites- 
J jt j+2 ; 
() 


It is apparent that the remaining terms have been reduced by 
factors equal to those obtained in the periodic case. 

The above result has been verified experimentally by process- 
ing a periodic signal containing only a strong bottom reflection 
(R=.55) and multiples (figure 47a). Removal of the first 
multiple component in the cepstrum results in 75% and 85% energy 
reductions of the second and third multiples respectively. 

These correspond exactly to the 50% and 67% amplitude reductions 
predicted by (11). The periodicity was then disturbed by 
average deviations of 1%, 5%, 10% and 20% of the original 
period. In each case the later multiple reductions coincided 
with (11). Figure 47 shows waveforms before and after process- 


ing for the exactly periodic signal and the case of 20% 


mPewlOdL Cl ty . 
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(a) 


Figure 47 (a) Sigral composed'of only periodic multiples of 
the bottom response. 
(b) Result of removing the first multiple cepstrum 
GOnER bu ELOn. 
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(c) 


(d) 


- Figure 47 (cont'd) (c) Signal of (a) with periodicity disturbed 
by average deviations of 20% (.2 sec). 
(d) Result of removing the first multiple 
cepstrum contribution. 
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3. Multiple-to-Signal Ratio 

The results obtained for aperiodic reverberation must be 
extended to signals containing reflectors as well as multiples 
if they are to be of practical use. This is not feasible 
analytically due to the nonlinearity of the homomorphic trans- 
formation and the reflector characteristics in the earth model 
used. Recall that, for this analysis, reflectors have the 


form t“e Pt 


rather than simply impulses. Each seismogram con- 
tains an additive combination of these reflectors with 
multiples. The logarithmic operation on the Fourier transform 
of this sum causes the cepstral contributions due to reflectors 
and multiples to be analytically indistinguishable. Hence, 
this phenomenon has been investigated empirically for various 
reflector sizes. It was found that cepstral properties of 
multiples are essentially preserved in the presence of non- 
impulsive reflectors although important effects were evicent 
in the cases tested. 

Percentage removal of the first multiple was found to be 
dependent on the width of the cepstral stopband. Figure 48 
illustrates this effect for seismograms with different reflector 
stengths and periodic multiples. Each signal requires a 
stopband of about 200 msec (41 samples) for complete removal 
Bamtme first multiple. As notch width 1S decreased the pertor— 


Mance becomes increasingly sensitive to MSR. For the smallest 


notch width shown (40 msec), performance varies almost 30% 
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Figure 48 Effect of cepstral stopband width on homomorphic 
dereverberation performance for signals of varying 
MSR. 
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with MSR. In the absence of reflectors the first multiple 

eould be completely removed by zeroing one sample of the cepstrun. 
It appears that the inclusion of reflectors has the effect of 
Sistributing the multiple energy in the quefrency domain, 
although, in all cases tested the energy remains concentrated 
Near the time domain multiple location. Figure 49 shows the 

0.5 second section of the cepstrum centered at the multiple 
location for each of the seismograms of figure 48. The cepstra 
are identical in form but the absolute cepstral energy increases 
with MSR. Three large peaks occur in each cepstrum between 

.94 and 1.0 second. This similarity in form suggests that 

Boal Bp ands should produce equivalent percentage reduction 
of multiples. We conclude, Mee ere from the experimental 
results that a greater portion of the multiple energy is 
Concentrated near 1.0 second in the higher MSR cepstra. This 
effect makes it difficult to select stopband limits by peak 
detection or visual inspection of the cepstrum. NOtch widths 
which yield the meee trade-off between multiple reduction and 
merlector distortion must be HeLa by trial-and-error 


Pomeparticular applications. 


4. Multiple Distortian 
Since homomorphic dereverberation is theoretically based 
on a Model of strictly impulsive multiples, it is important to 


observe the performance of this technique in the more realistic 
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Figure 49 0.5 second section of cepstra at multiple location. 
(a) MSR=20 (b) MSR=1 (c) MSR=.06 
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case of an extended bottom interaction mechanism. We employ 
the same bottom response functional form discussed earlier 

in this chapter. The seismogram of figure 50a contains 
distorted multiples due to a bottom response of .045 seconds 
duration. Reflectors are present at 2.7 and 3.5 seconds. 
Figure 50b shows the region of the cepstrum centered at the 
first multiple location. The peaks are much broader than those 


of figure 49. The dominant energy is concentrated near the 


multiple location as before but it now appears to be more 


distributed. 

Table 1 summarizes the effects of varying the stopband 
width and location in the cepstrum of figure 50b. First and 
Phird multiple energy removed varies widely while the second 
multiple is not significantly affected by the passband dimensions. 
In some cases, extending the stopband decreases the amount of 
multiple energy removed. Zeroing the .94-.98 second region 
results in greatly increased reduction of the third multiple 
(74-783) but first multiple reduction is degraded by 10-12%. 
This behavior, which has been observed in several cases, is 
not predicted by the theory and is thought to be a eonoucae one! 
effect. The high amplitude oscillations which are dominant in 
the left side (.75-1.0 second) of figure 50b, or the low 
frequency "drift" which is apparent throughout the section may 


be computational noise which contributes to this phenomenon. 


Several observations can be made in this case. Although 
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(a) Seismogram with distorted multiples 
(b) 0.5 second section of cepstrum Semis Tete at 
first Multi plLesvoca elon, 
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Stopband Limits ENERGY REMOVED 
(sec) 
lst MULT 2nd MULT _3rd MUL Ist REPL aezne Rei. 


pool .04 isis 


-99-1.04 | 39 | -.12 -.05 | 
a 
| 





1 
© 
oO’ 


1.02-1.06 45 [| .50 | -.13 





TABLE lL 


EFFECTS OF STOPBAND LIMITS ON PERFORMANCE 
FOR THE SIGNAL OF FIGURE 50. 
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performance is somewhat erratic, a stopband of 40-80 msec 
which straddles the first nee ioe onset time yields significant 
Multiple reduction. The peak performance in this case is about 
10% below that achieved with undistorted multiples but the 
sensitivity to minor stopband variations is less dramatic. 
Reflector distortion does not increase significantly due to the 
presence of distorted multiples. 

Figure 51 shows processed seismograms resulting from the 
application of various stopbands to the cepstrum of figure 50. 

Finally, we note that the effects of aperiodicity have 
not been discussed in the distorted multiple case. Due 28 
limitations of the earth model it was not possible to investigate 
these effects. Thusfar, however, we have relaxed the theoretical 
assumptions regarding periodicity and coherence individually, 
and found that homomorphic dereverberation is not critically 
sensitive in either case. We surmise that the peered Or 
these effects would not be catastrophic to performance. Further 
investigation is merited since this topic has an important 
bearing on the effectiveness of the homomorphic technigue in 
shallow water dereverberation where later multiples are 


Significant. 


5. Water Travel Time Estimate 
In practice, passband location must be determined by 


estimation of water column travel time. It is apparent from 
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Figure 51 Processed results of the seismogram in figure 50a 
for three different cepstral stopbands. (a) 1.02-1.06 sec 
(b) 1.0-1.04 sec (c) .98-1.06 sec 
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the preceding discussion that the importance of accurate 
bottom estimation depends on the seismic environment. In 
particular, we have seen that MSR and bottom characteristics seen 
memaetect cepstral energy distribution. Proximity of multiples 
tO important reflectors also dictates resolution constraints 
on travel time estimation. For the cases analyzed here the 
required travel time estimation accuracy would be 20-40 msec 
Since, as determined in the foregoing discussion, a stopband 
of 40-80 msec is generally required for effective dereverberation. 
Travel time estimation error of more than half the stopband 
width causes the major multiple contributions to be excluded 
from the stopband. This degree of accuracy can easily be 
attained with existing bottom eee ets algorithms. 

The reflector/multiple resolution implied by the required 
stopband widths is also about 40-80 msec for the data tested. 
A discussion and several examples of reflector/multiple resolu- 


tion are included in the following section. 


6 Reflector/Multiple Overlap 

The following six figures illustrate homomorphic dereverbera- 
Bion Of Signals in which reflectors and multiples are closely 
Situated. 

Figure 52a shows a seismogram with reflectors at 1.6, 2.02 
and 3.0 seconds. The first and second multiples directly inter- 


fere with reflectors and the third multiple is barely visible 
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(a) 


(b) 


Figure 52 (a) Unprocessed seismogram 
(b) Result of cepstral notch filtering; .96-1.02 sec 
stopband 
(c) .98-1.06 sec stopband 
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at 4.0 seconds. The result of applying a cepstral stopband at 
.96-1.02 seconds is seen in figure 52b. Most of the energy 
near mre) seconds has been TeMmoVea {oUer ay smal ipe rec ielemS 
Still visible. The signal of figure 52c results from a stop- 
mama Of .98-1.06 seconds which spans the time domain reflector 
location as well as the multiple location. Residual energy 
1S still present after this processing which indicates that 
some of the reflector and first multiple cepstrum contributions 
are distributed beyond the immediate vicinity of the time 
domain locations. 

The seismogram of figure 53a contains reflectors at 1.6, 
1.95 and 3.0 seconds. Each reflector is clearly evident after 
Processing (.98-1.06 second stopband) and first multiple energy 
has been greatly reduced as shown in figure 53b. Some of the 
removed energy may be due to the 1.95 second reflector; however, 
the first multiple in this example is very coherent and its 
energy is more likely to be localized near 1.0 second in the 
cepstrum. Extension of the stopband closer to the reflector 
location (.96-1.06 second) causes complete removal of the 
reflector as seen in Figure 52C.9) Visible aedte eon Ou jae 
Second multiple at 3.0 seconds and an internal multiple at 2.7 
seconds is apparent in figure 52b. 

The seismogram of figure 54a has the same reflector place- 
ment as that of figure 53a but the multiples in this case are 


Smaller and less coherent. Application of a cepstral stopband 








Figure 53 
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(a) Unprocessed seismogram (hb) -Result of cepstral 
notch filtering; stopband .98-1.06 sec. (c) .96-1.06 
sec stopband. 
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Figure 54 (a) Unprocessed seismocram . . 
: (pb) Result of cepstral nNoOtechetititeuwma Stoppard 
1.0-1.06 sec. 
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at 1.00-1.06 second yields (figure 54b) visible reduction 
of all three multiples and well resolved reflections at the 
proper locations. : 

Figure 55a illustrates a seismogram with considerable 
multiple distortion and reflectors at 1.5, 2.1 and 3.0 seconds. 
A wide cepstral stopband, .98-1.08 second, has a very minor 
effect on the first multiple as seen in figure 55D, EXtens ton 
of the stopband to 1.12 seconds causes removal of virtually 
all signal energy in the region. The relatively wide (.1 second) 
reflector/multiple separation cannot be effectively exploited 
in this case because the very incoherent multiple in a low MSR 
signal requires a large stopband for effective removal (see 
figure 48). 

The result of longpass filterina the cepstrum of the signal 
faeriagure 55a is shown in figure 56. All cepstral contributions 
prior to 1.02 seconds have been set to zero. The bottom and 
later reflectors are clearly visible but the 1.5 second reflector 
has been removed. This illustrates one drawback to longpass 
dereverberation in deeper water. This effect may be acceptable 
in some situations, however. For instance, good resolution of 
closely spaced, deep reflectors may be obtained by longpass 
filtering whereas the earlier reflectors are frequently obvious 
before processing. 


Figure 57 shows longpass results foOrva sigeah wlth Impulsive 


multiples and very sharp reflectors at 2.2 and 3.0 seconds. 
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Figure 55 


(a) 
(b) 


(c) 


Unprocessed seismogram. 

Result of cepstral notch filtering; stopband 
2G en ome ee 

-98-1.12 sec 
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(a) 


Figure 56 (a) Unprocessed seismogram 
(b) Result of longpvass filtering cepstrum; cutoff 
gquefrency 1.02 sec. 
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(a) 


Figure 57 (a) Unprocessed seismogram 
(b) Result of longpass filtering; cutoff quefrency 


1.01 sec. 
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The stopband limit is 1.01 seconds in this case. Multiple 
removal is complete and both reflectors are clearly visible. 
Longpass filtering was found to be most effective for signals 
of this type. Very sharp reflectors are more clearly visible 
than more distorted reflectors of comparable energy, especially 
in noisy signals. lLongpass filtering of noisy signals is 
discussed later in this chapter. 

The low frequency noise present in the reflector region 
Preergure 57 is typical of that observed in several longpass 
results. The cepstra of figures 49 and 50 contain similar 
components. No attempt was made in this analysis to remove 
Bats type of noise. The reason for its presence has not been 


determined. 


7. Additive White Noise and Filtering 

The effects of additive noise are not addressed in the 
formulation of the homomorphic system for deconvolution since 
there are no avvarent cepstral properties of noise which can be 
exploited for its removal. Linear filtering is a more suitable 
way of reducing additive noise in individual signals. Thais sean 
be Bee rorined orior to, or in conjunction with, homomorphic 
deconvolution. The effects of this combined processing in the 
special case of Gaussian white noise have been investigated and 
are discussed here. The data presented represent a relatively 


small number of experiments performed with synthetic data and 
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noise generated as previously discussed. These results cannot 
be generalized categorically due to the relatively narrow 
scope of the experiments and the lack of precise theoretical 
characterization of noise properties under homomorphic trans- 
formation. The results do indicate performance trends and 
computational effects due to additive noise and filtering. 

Since noisy Signals usually undergo linear filtering prior 
to dereverberation, we begin by discussing an important 
computational issue which arises when signals are bandpass 
filtered. Such filtering introduces spectral regions (stopbands) 
containing little or no spectral energy. Recall that the homo- 
Pemehic meet ormetion involves computing the logarithm of the 
Fourier transform of the signal. Simee the logarithm of zero 
is not defined, this computation is not possible in frequency 
bands where the Fourier transform is zero. In the case of 
lowpass filtered signals this problem can frequently be overcome 
by resampling after filtering. Figures 58 a, b and c illustrate 
this process. If the sampling frequency is decreased to the 
Nyquist rate implied by the filter cutoff frequency, the 
resulting ueerare Fourier transform will: include only the 
frequency components in the passband region. Figure 58d shows 
a distribution of spectral energy which is not readily amenable 
to elimination of the zero region. In this case the baseband 
region is zero so that resampling would not be effective. 
Investigation of this problem is currently in progress 


(Tribolet [11]). 
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Figure 58 (a) Unfiltered spectrum sampled at the Nyquist rate. 
(b) Spectrum of (a) after lowpass filtering at WwW. 
without resampling. 
{c) Filtered spectrum after reampling. 
(d)} Spectrum with zeros 1m Ene escbanaeregton 
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it was found in this analysis that computing the cepstrum 
of an over-sampled Signal is an unstable procedure which, at best, 
requires extensiye computation time. The low energy regions 
Pemene Spectrum lead to spurious phase derivative values. 
Hence, the phase unwrapping algorithm proceeds at small step 
sizes, requiring computation of many intermediate values of the 
discrete Fourier transform. In many cases the computer algorithm 
Sould not produce an unwrapped phase which was in acceptable 
agreement with the principal value. Cepstra of over-sampled 
Signals which were successfully computed, however, generally 
yielded dereverberation results comparable to those of properly 
Sampled signals (see table 2). 

Addition of white noise was found to have a computational 
effect similar to that described above; phase unwrapping time 
is significantly increased. Heavy weighting (w =.98-.99) 
was found to reduce computation time considerably for noisy 
signals, although some signals require small phase integration 
step sizes in isolated sections even when substantial weighting 
is applied. Recall from Chapter II that there is no assurance 
that the log-spectrum is adequately sampled, even after lowpass 
Mmiikeering of the signal. Hence, the integration Gf the phase 
derivative cannot be expected to prcceed quickly in all cases. 

Smoothing of the phase derivative was attempted to 
compensate for the effects of noise. A three-point moving 
average was applied prior to integration. Themsesultingsinvers 


cepstra bore no resemblance to the Original seismograms. 
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TABLE 2 


EFFECT OF RESAMPLING ON MULTIPLE REMOVAL FOR A 
NOISELESS SEISMOGRAM LOWPASS FILTERED AT 50 and 

20 Hz. DEREVERBERATION WAS ACCOMPLISHED BY APPLY- 
ING AN 80 MSEC CEPSTRAL STOPBAND AT THE FIRST 
MULTIPLE LOCATION. 
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Figure 59 summarizes the observed effects of noise on 
first multiple removal by homomorphic processing. Each 
Seismogram was lowpass filtered at 50 Hz and resampled by a 
factor Of four prior to multiple removal. Percentage of 
multiple energy removed shows a consistent decrease with increas- 
ing noise level. The rates of decrease and amounts of energy 
removed are seen to vary widely from one signal to another. 
Examples of noisy signals before and after processing are shown 
imeragures 60 and 61. Significant reduction of the first multiple 
is evident in both examples. In the first case, the noise level 
is moderate and multiple reduction has a marked effect ef wileuen 
quality of the signal. The noise level in figure 61 is consider- 
ably higher, resulting in marginal improvement due to dereverbera- 
j10On. 

Second and third multiple energy removed was found to 
decrease generally with decreasing noise also, and in some 
high noise cases the later multiples were actually enhanced. 


Data for a typical signal are tabulated below. 


STANDARD 


DEVIATION ENERGY REMOVED 


OP NOLS 2nd MULT 3a" Muna 





TABLE 3 
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Figure 59 Effect of added noise on homomorphic dereverberation. 
All signals are lowpass filtered at 50 yg and 
resampled. 
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i = Hz. 
(a) Seismogram with very high noise, F_=50 
Figure 61 (pb) Result Of cepstrommuoren filtering’ (.80-1.04 sec). 
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Reference to figure 61 reveals that for the higher noise 
amplitudes the later multiples are buried so that their removal 
is not important. 

Reflector distortion, summarized for a CY pled wecase win 


table 4, was found to be generally more severe than in noise- 


less Signals but rarely more than 25%. 


STANDARD 


DEVIATION REFLECTOR ENERGY REMOVED 
Or NOISE l 2 2 


eae [ee 


50 7 =e led OS 
100 14 =,/2 Seno 
TABLE 4 


EFFECT OF NOISE ON REFLECTOR DISTORTION 
FOR A TYPICAL SEISMOGRAM WITH 3 REFLECTORS, 
FILTERED AT 50 Hz. 


The above results indicate that the addition of noise 
leads to decreasing performance with respect to all three 
criteria. The behavior is somewhat spurious and frequently 


exhibits wide variation from signal to signal. 
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several lowpass filter bandwidths (2055 30, 490; JOP and 
90 Hz) were applied to the signals evaluated in figure 59. 

The results are shown in figure 62. In each case, first multiple 
removal is least effective when the signals are prefiltered at 

20 Hz. Figures 63 and 64 illustrate this behavior. The same 
seismogram is shown in figures 63a and 64a, lowpass filtered 

at 530 and 20 Hz respectively. The difference in appearance is 
dramatic. Resolution of the reflectors on either side of the 
fest multiple (at 1.55 and 2.4 sec) is greatly improved by 
dereverberation in the latter case while figure 63 shows very 
little visual improvement. 

The behavior illustrated in these figures cannot be fully 
explained on the basis of the Aiea available. The signals 
tested have dissimilar reflector locations and varying amounts 
of multiple distortion which implies that the similar performance 
dips are not due to similarities in signal configuration. The 
observed behavior may be due to the properties of the source 
Signature (which is common to all three seismograms) or the 
characteristics of the recursive third crder Butterworth 
filter employed. More data using different filter routines 
and a wide range of signal characteristics is needed to 
completely characterize this Hehe The results obtained 
here suggest that a filter bandwidth very close to the bandwidth 
of the noiseless signal leads to the best dereverberation 


performance. The same result was obtained for the TDL algorithm. 
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Figure 62 Effects of lowpass filtering signals with added 
white noise prior to homomorphic dereverberation. 
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Figure 63 (a) Noisy seismogram lowpass filtered at 50 Hz. and 
resampled at 51 Hz. 
(b) Result of ceostral notch filtering; stopband 
.8-1.04 sec. 
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Figure 64 Noisy seismogram lowpass filtered at 20 Hz and . 
resampled at 51 Hz. (b) Result of cepstral filtering; 
stopband .8-1.04 sec. 
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Additive noise was found to have a particularly adverse 
effect on longpass filtered signals. Noise, which may be 
relatively low in the received signal, tends to be amplified 
heen respect to reflectors in the process of longpass filtering, 
especially in the later portion of the signal. Although long- 
pass filtering removes a large part of the noise energy with 
the source, the overall effect is usually a decrease in SNR. 
The advantage of longpass filtering, which includes source 
deconvolution, is that greater resolution of close reflectors 
can be achieved. It was found that this approach is worthwhile 
in low noise signals but not effective when the white noise 
level is significant with respect to reflector amplitudes. 

Proper resampling after prefiltering was found to be very 
important for successful longpass filtering. Figure 65a 
illustrates the effect of longpass filtering the cepstrum of a 
noisy signal which was first lowpass filtered at 20 Hz but 
not resampled. The filtered signal contains relatively low 
noise and the sampling frequency is 205 Hz. The processed 
result of figure 65b is useless due to the high sampling rate. 
Reduction of the sampling frequency to 102.5 Hz leads to the 
processed signal of figure 66a. The 3.5 second reflector is 
clear but the high background noise almost obscures the 2.7 
second reflector. Multiple removal is complete. Resampling 
to 51 Hz, which is approximately the Nyquist rate in this case, 


leads to some improvement (figure 66b) but the SNR is much lower 
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Figure 65 (a) Noisy seismogram lowpass filtered at 20 Hz but 
not resampled. 
(b) Result of longpass filtering the cepstrum of (a). 
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Figure 66 Longvass processing results (a) for the signal of 
figure 65a, resampled at 102.5 Hz before vrocessinga. 
(b) for the signal of figure 65a, resampled at 51 Hz 
before orocessing (c) for a signal identical to 
figure 65a with one half Ehe wn teemmoree fever 
resampled at 51 Hz. 





a 3= 


than before dereverberation (figure 65a), Figure 66c results 
from processing identical to that of figure 66b, on the same 
Signal, with the noise amplitude halved. Both reflectors are 
Clear and the multiples have been removed but the background 
noise is much higher even than in the signal of figure 65a. 
Thus, longpass filtering of noisy signals is seen to involve 

a trade-off between effective dereverberation and decrease in 
SNR. For signals with moderate to heavy noise the reduction in 
SNR was found to be unacceptable in the cases tested. 

Although this subject was not extensively investigated 
there 1S some indication that lowpass filtering can be employed 
to improve the results of subsequent longpass processing. 
Figures 67 and 68 illustrate the longpass processing of a 
noisy signal after lowpass prefiltering at (a) 70 Hz, (b) 50 Hz 
and (c) 30 Hz. In each of the signals in figure 68 the multiples 
have been removed but the reflector resolution improves 
considerably from (a) to (c). These results are reasonable 
in that improvement of performance coincides with increasing 
rejection of out-of-band noise; however, the filter and signal 
characteristics must be studied more closely to explain this 


behavior accurately. 


D. Comparative Examples of Processing Results 


The foregoing results and discussion illustrate the 


performance of the TDL and homomorphic dereverberation algorithms 
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Figure 67 Noisy seismogram lowpass filtered at three different 
frequencies (a) 70 Hz (b) 50 Hz (c) 30 Hz. Each has 
been resampled at 51 Hz which is the approximate 
Nyquist rate for the signal without noise. 
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Figure 68 Resuits of longpass filtering the seismograms of 
figulle 65 a, b and c, Yrespecci ele 
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for a variety of processing conditions. Several trends in 
performance are apparent. Before proceeding to a summary 

and A gemgetion of the relative strengths and weaknesses, we 
present several examples which allow direct visual comparison 
of the techniques. Each of the following figures contains 
(a) an unprocessed seismogram and the two processed results 
Obtained by (b) homomorphic and (c) TDL filtering. 

Figure 69 illustrates a noiseless seismogram with multiple/ 
reflector overlap at both 2 and 3 seconds. Both algorithms 
leave asmall amount of energy at the first location and effect 
Only slight reduction at the second. In general, both methods 
were found to eradicate reflectors which are extremely close 
to the first multiple and retain signal components which closely 
coincide with later multiples. 

In figure 70 the multiple onset occurs .1 second before 
the reflector and the MSR is considerably lower than in the 
previous figure. In this case the separation is great enough 
that both methods retain a significant amount of signal energy 
near 2.1 seconds. The homomorphic result is considerably 
sharper although very little energy has been removed. Multiple/ 
reflector separation of .1 second was found to be the approxi- 
mate resolution limit of both techniques when reflector onset 
is later than multiple onset and travel time is estimated 


accurately. 


Figure 71 shows a reflector at 1.95 seconds, slightly 
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Figure 69 (a) Unprocessed seismogram with reflectors at 1.5, 2.02 


and 3.0 seconds. 
(b) Result of homomorphic processing. 
(c) Result of TDL process taee 
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Figure 70 (a) Unprocessed selismogram with reflectors at 125.0 222 
anGaone SEC. 
(b) Result of homomorphic processing 
(c) Result of TDL processing: 
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Unprocessed seismogram with reflectors at 1.6, 1.95 
Gone) seer 

Result of homomorphic vrocessing. 

Result of TDL processing. 
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before the multiple. Effective dereverberation is accomplished 
in both b and c. As in the previous figure the reflector is 
better resolved by homomorphic processing, Resolution of both 
techniques was generally observed to be slightly better when 
reflector onset is earlier than multiple onset, provided travel 
time is estimated accurately. In such cases of severe inter- 
ference the homomorphic filtering generally yields more distinct 
reflections. A further example of this behavior is shown in 
figure 72. 

The following three figures illustrate dereverberation of 
noisy Signals. Each seismogram has been lowpass filtered at 
30 Hz and the homomorphic outputs as shown have been resampled 
fMemolL HZ. in this first case (figure 73) the performance of 
both methods is comparable. The first multiple is almost 
completely removed while other regions of the signal are not 
visibly affected. Figure 74 also shows comparable multiple 
removal, however, the resolution of the reflectors at 2.6 and 
3.4 seconds is somewhat better after TDL filtering. The 
homomorphic algorithm was found to produce higher random noise 
Spikes than the TDL filter. This effect is present in figure 
74 and again in figure 75. In both examples the homomorphic 
method achieves slightly better multiple removal but the 


overall noise level in the result appears higher. 
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Unnrocessec seismogram with reflectors at 1.6, 
1295 ane 2.0) secs 

Result of homomorphic processing. 

Result of TDL processing. 
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Figure 73 (a) Noisy seismogram with reflectors at 2.7 and 


3.5 sec, and lowpass filtered at 30 Hz. 
(b) Result of homomorphic processing. 
(c) Result of TDL processing. 
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Figure 74 (a) Noisy seismogram with reflectors at 1.55, 2.4 and 
3.25 sec, lowpass filtered at 30 Hz. 
(b) Result of homomorphic processing. 
(c) Result of TDL processing. 
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Figure 75 (a) Noisy seismogram with reflectors at 2.6 and 3.4 


sec, lowoass filtérea@ar siesta 
(b) Result of homomorphic processing. 
(c) Result of TDL processing. 
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CHAPTER V 


DISCUSSION AND SUMMARY 


In the preceding chapters we have (1) examined the 
theoretical structures of the TDL and homomorphic dereverbera- 
tion techniques, (2) established a comprehensive set of perfor- 
mance criteria, and (3) presented the results of applying both 
methods to synthetic data. Our approach has been essentially 
that of perturbation analysis. Through variation of environ- 
mental and signal processing parameters we have observed 
performance trends due to deviations from the ideal theoretical 
models upon which the methods are based. In a more practical 
sense the parameter variations simulate a range of seismic 
processing conditions. Since there are many different environ- 
ments in which these algorithms may be aaa. we have not 
emphasized the absolute performance figures obtained from the 
Simple, synthetic data utilized here. Rather, we have tried 
to present a behavior profile of both algorithms which indicates 
the basic trends and sensitivities with respect to a number 
of parameters, interpreted in terms of their theoretical 
structures and assumptions. This approach is intended to give 
Aemore general indication of the dereverbewartonepercomance 
to be expected in different situations. We conclude with a 


summary and discussion of the comparative results presented in 


Chapter IV. 





| 


| 


= 56 = 

In low noise, deep water seismograms the methods have 
been found to be comparable for reducing dominant first 
multiples (by 75-95% in most cases). The TDL filter requires 
accurate Signal statistics including an estimate of the multiple- 
reflector crosscorrelation function, which must be approximated. 
This leads to degraded performance in shallow water situations 
where Significant reflector energy is within the crosscorrelation 


window. The homomorphic method requires no statistical character- 


ization of the signal and thus has no similar performance 


degradation in shallow water; however, the three-stage cepstral 
transformation requires extensive computation which may be an 
important limitation for at-sea processing systems. (This 
issue will be discussed in more detail later.) 

Although the effects of aperiodicity Coula not be theroughiily 
evaluated experimentally the derived result expressed in 
equation (11) suggests that the homomorphic algorithm has the 
potential to reduce later, aperiodic multiples. The combina- 
tion of such processing with source deconvolution by longpass 
filtering the cepstra of shallow water seismograms appears to 
be the most promising application of homomorphic dereverbera- 
eon . Mere extensive practical evaluation is needed in this 
area. 

Closely spaced, aperiodic multipllesmaestue, ste conerenee 


of the approximated crosscorrelation function at shifts near 


the two-way travel time which, again, limits the effectiveness 





See 


of the TDL filter in shallow water signals. 

Increasing white noise level causes monotonically decreas- 
ing performance in boca Ehage as illsustrated in figures 38 
and 57. The TDL performance fell off more slowly with noise 
when both techniques were applied to.similar signals. Bandpass 
filtering leads to a consistently higher percentage of multiple 
reduction by TDL processing of noisy signals. Filter effects 
On homomorphic processing are more complex. The cases evaluated 
indicate that multiple energy removed is not a monotonic func- 
tion of filter cutoff frequency. Considerably more data are 
required to determine the precise effects of filter bandwidth 
and phase characteristics. The resampling which was found to 
be helpful after bandpass filtering may have a detrimental 
effect on visual record quality, so that interpolation may be 
desirable in some cases. 

Reflector distortion does not appear to be a problem for 
either technique except in cases where overlap is severe. In 
most cases tested less than 10% of the reflector energy was 
removed. Close proximity of reflectors to the first multiple 
frequently leads to severe distortion by. both methods due to 
the resulting bias of the crosscorrelation function and the 
lack of sufficient cepstral separation. Reflectors close to 
later multiples are usually well preserved. This behavior 
suggests two applications of the homomorphic algorithm. First, 


when regions of geological interest occur after the onset of 





=o 


the first multiple, a wide cepstral stopband can be employed at 
the first multiple location to reduce later multiples without 
distorting reflectors. In very shallow water the stopband may 
be extended to include more than one multiple. A second 
possibility for avoiding reflector distortion is the use of 
weighting coefficients greater than 1.0 to exploit the properties 
of mixed phase sequences. The object of this weighting is to 
make the reflector train mixed phase while keeping the multiple 
sequence minimum phase, This appears to be feasible in many 
Situations since the z-plane zero of the multiple sequence is 
usually well inside the unit circle. Moving some of the reflector 
train zeros outside of the unit circle (i.e., making it mixed 
phase) will, in general, cause some of the cepstral energy due 
to the reflectors to occupy the negative quefrency region. 
Even if the reflector train has a eee phase component 
before weighting the same effect can be expected. Thus, the 
amount of reflector energy near the first multiple location 
may be reduced. Although there is no guarantee that the 
resulting notch filtered cepstrum will transform to a seismogram 
with less distortion, this technique appears to be worthy of 
investigation. 

We recall one other reflector distortion effect which was 
seen in Chapter IV. We saw in figure 55 that dereverberation 
by longpass filtering completely removes reflectors eccurr im 


prior to multiple onset. It was noted that this effect may be 
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acceptable if it leads to better resolution of smaller, deep 
reflectors. 

In terms of our third criterion, visual improvement of the 
Signal, both techniques were seen to have advantages and 
disadvantages. Homomorphic processing usually resulted in 
better resolution of interfering signal components in the low 
noise signals processed. Increasing noise, however, was seen 
to cause homomorphic results more prone to random noise spikes 
which degrade the interpretability of the record. The results 
of TDL filtering had somewhat better visual resolution in the 
noisier signals processed. As noted previously the visual 
advantage of the TDL in this case is partially due to the 
noisier appearance of the rosamened Signals produced by 
homomorphic processing. 

Longpass filtering was Seen to provide the most effective 
dereverberation, the best reflector resolution and the best 
overall visual quality in ideal cases. Unfortunately it 
degenerates quickly with noise and could not be successfully 
applied to very noisy signals or signals with important geological 
regions preceding the multiple onset. Further research Baa 
experience may well lead to more extensive applicability of 
this technique. | 

The relative simplicity of the TDL algorithm makes it 
much more desirable from a computational standpoint. TDL 
dereverberation of a 1024 point signal can be accomplished in 


3 seconds or less for the operator lengths typically required. 
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The three major computational steps are the correlation 
operations, solution of the estimator equations and convolu- 
eon Of me operator tee the signal. 

The homomorphic computations include weighting, four FFT's 
computation of the complex legarithm, phase unwrapping, linear 
filtering, complex exponentiation and unweighting. This 
algorithm can be expected to take 20 seconde more for a 
1024 point sequence on a small processing computer. In this 
analysis the phase unwrapping computation took over one minute 
for some noisy signals. These figures are highly dependent 
on hardware available and programming efficiency but, in 
general, homomorphic dereverberation is several times slower 
than the TDL algorithm. Special purpose hardware could be 
used to reduce homomorphic computation time significantly, 
but the method has not been implemented for processing on a 
large scale thusfar. 

Storage requirements for the homomorphic algorithm vary 
with the FFT routine used, method of cepstrum computation and 
cepstrum length. The program used for this analysis requires 
about 12 * N bytes of core and 4 * N bytes of disc storage, 
where N is the cepstrum length. N was twice the signal length 
de the cepstra computed. Shorter cCepstrunm@tengens, as 
determined by trial-and-error, may produce results with 


acceptably low aliasing in many cases. The TDL dereverberation 


program used requires about 6 * L bytes of core, where L is 
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the data sequence length. No disc storage is required in 
this computation. These sconaae requirements apply to a float- 
ing point processing scheme on a machine (HP-2100) which uses 
four byte floating point words. 

In conclusion, we make some general observations concern- 
ing the results of this analysis. 

The TDL dereverberation scheme is a simple and efficient 
technique which has demonstrated effectiveness in removing 
deep water multiples. The analytical structure is well under- 
stood and its performance characteristics have been explained 
here in terms of that structure. Further refinements in its 
implementation may be possible but its potential is essentially 
clear at enis Ole ak quae 

Homomorphic dereverberation 1s complex, relatively untested 
and requires extensive computation. It has been shown here 
to be effective on synthetic data. It appears to be particularly 
promising for shallow water dereverberation. The complexity 
of the method leads to a number of possible analytic and 
computational techniques which can be utilized in its applica- 
tion. Certainly, its full potential has not yet been determined 


and further investigation is warranted. 
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APPENDIX A 
COMPUTATION OF THE PHASE DERIVATIVE OF THE Z-TRANSFORM 

When computing the complex cepstrum of a sequence, x(n), 
it is necessary to determine the unique, continuous phase of 
X(z). One way of obtaining the continuous phase is to first 
compute its derivative and then integrate eee The 
computation of the phase derivative from x(n) is discussed 
in detail here. 


We begin with the z-transform of x(n), 


oo 


X(z) =) x(n) 2 exe (ee ence 


=—0co 


Taking the complex natural logarithm 
log X(z) = X(zjy"= logul eG ree ce) 
We see that the phase of X(z) is equal to the imaginary 
Bart Of Ets natural logarithm. The wdenivaetye= son X,(Z) can be 


expressed in terms of easily computable quantities. 


a2) ee cate ee ri 
dz dz X(z) 


where the prime indicates differentiation with respect to 2a. 


Expanding (1) and solving for X, (2), 
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X# (2) + 4X} (2) 


X'(z) + 4X. (2) 
= z Xp (z) + JX, (2) 


A Xi (2) = (2) 


>< 

vx 

N 
ll 


+ 3X4 (2) 
X, (2) + 3X, (2) 


pePDarating the RHS into real and amagimany, parts) 


e Ke (2) he (2) te eee 
Ki (2) = R I I R 


eat r x (2) 


(X,(z) Xi (z) + X(z) Xi (z)) 


oa 
y= 


7 (z) - 


Ke (2) + xe (z) 


The real part yields an expression for the phase derivative, 


wy) Xyl2) Xp (2) ~ Xz (2) x, (2) 
ee + Soe 
R I 


Since the z-transform is actually evaluated on the unit circle 


using the discrete Fourier transform (DFT) we set z = eJ%. 


a x (eI) ee - X (eo x' (eI) 
Stel) = 2 eee (3) 


x2 (24) + x2 te!) 


2 


Derivatives with respect toe : may be replaced by = since 
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e s f e 
dx(e7")| _ 5 giv jax(e?) 
du deJ™ 


and we thus have a common factor of j e/” in each term of (3). 

Hence, we need only compute the real and imaginary parts 
of x (eI) and = (x (eI) ) and combine them as indicated in (3). 
The derivative of x (e)%) is easily obtained from the sequence 


n x(n) as follows: 


oo Le Jw 
os) Se ee 
Re [X_(e7") ] [ ~x1(eJ") 
Im[x, (e2%)] = x ee 


The required computation in terms of the DFT is then 
a4 -(x,(k) X p(k) + X,(k) Xy (e)) 
X_(k) = SC 
; x2 (k) + x2 (k) 
Rea ue 
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